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_^ ! Abstract 
in 

, The general relativistic corrections in the equations of motion and associ- 

. ated energy of a binary system of point-like masses are derived at the third 

. post-Newtonian (3PN) order. The derivation is based on a post-Newtonian 

O . expansion of the metric in harmonic coordinates at the 3PN approximation. 

I , The metric is parametrized by appropriate non-linear potentials, which are 

$—1 . 

evaluated in the case of two point-particles using a Lorentzian version of an 
Hadamard regularization which has been defined in previous works. Distri- 
^ , butional forms and distributional derivatives constructed from this regular- 

ization are employed systematically. The equations of motion of the particles 
are geodesic-like with respect to the regularized metric. Crucial contribu- 
tions to the acceleration are associated with the non-distributivity of the 
Hadamard regularization and the violation of the Leibniz rule by the distri- 
butional derivative. The final equations of motion at the 3PN order are in- 
variant under global Lorentz transformations, and admit a conserved energy 
(neglecting the radiation reaction force at the 2.5PN order). However, they 
are not fully determined, as they depend on one arbitrary constant, which 
reflects probably a physical incompleteness of the point-mass regularization. 
The results of this paper should be useful when comparing theory to the 
observations of gravitational waves from binary systems in future detectors 
VIRGO and LIGO. 
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I. INTRODUCTION 



The present work is a contribution to the problem of the dynamics of two compact objects 
at the so-called third post-Newtonian (3PN) approximation of general relativity. By 3PN 
we mean the relativistic corrections in the binary's equations of motion corresponding to 
the order 1/c^ relatively to the Newtonian acceleration, when the speed of light c tends to 
infinity. Why studying the equations of motion to such a frightful post-Newtonian order? 
A side reason is the strange beauty of the post-Newtonian expansion, which becomes quite 
intricate at the 3PN order, where it requires some interesting mathematical methods. The 
main reason, however, is that inspiralling compact binaries, namely systems of two neutron 
stars or black holes (or one of each) moving on a relativistic orbit prior to their final merger, 
should be routinely observed by the gravitational- wave detectors LIGO, VIRGO and their 
fellows. Several analysis show that the post-Newtonian templates required for the detection 
and parameter extraction of inspiralling compact binaries should include the relativistic 
corrections in the binary's orbital phase at approximately the level of the 3PN order 

Lorentz and Droste |^ were the first to obtain the correct equations of motion of two 
non-spinning particles at the IPN approximation (see for reviews). An important 

work by Einstein, Infeld and Hoffmann [p!T|-[l3| showed that the IPN acceleration can in fact 
be deduced from the vacuum gravitational field outside the masses. This result is interesting 
because, in their approach, the bodies are allowed to carry a strong internal gravity. Unfor- 
tunately, the computation of the surface integrals surrounding the masses is very difficult 
even at the IPN order (see [H] for a recent derivation of the Einstein- Infeld- Hoffmann equa- 
tions). The same equations were also obtained by Fock and followers [Jl5|-0 for the motion 
of the centers of mass of bodies with finite size. The next approximation, 2PN, has been 
tackled by Otha, Okamura, Kimura and Hilda |T^-^ with a direct post-Newtonian compu- 
tation of the Hamiltonian of point-particles; however the first complete 2-particle case in 
their framework is only given by Damour and Schafer and the fully explicit 3-particle 
case is due to Schafer [^. Up to the 2PN level, the equations of motion are conservative 
(existence of ten conserved quantities, including a conserved energy). The non- conservative 
effect, which is associated to the radiation reaction force, arises at the 2.5PN order. The first 
correct equations of motion of two masses at the 2.5PN order were obtained by Damour, 
Deruelle and collaborators [p^-p6| in harmonic coordinates. These equations are applicable 
to systems of strongly self-gravitating bodies such as neutron stars (see Damour P,pi]| for 
the proof). Moreover, Kopejkin and Grishchuk and Kopejkin |2^ obtained the same 
equations in the case of weakly self-gravitating extended bodies. The corresponding result 
at 2.5PN order was also derived by Schafer [^,0 using the ADM Hamiltonian approach. 
Later, the harmonic-coordinates equations of motion were re-computed by Blanchet, Faye 
and Ponsot |^ following a direct post-Newtonian iteration of the field equations. Some 
of the latter derivations |23,29,31] opt for a formal description of the compact objects by 
point-particles. There is a nice agreement between all these different methods at the 2.5PN 
order. In addition, the complete 2.5PN gravitational field generated by point-particles in 
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harmonic coordinates was derived in IST 



At the 3PN order, the equations of motion have been obtained using a Hamiltonian and 
formal delta functions by Jaranowski and Schafer P^ , PB| in the center-of-mass frame, and 
by Damour, Jaranowski and Schafer |^| in an arbitrary frame. These authors found an 
irreducible ambiguity linked probably with an incompleteness in the regularization of the 
infinite self- field of the particles. In this paper, following the method initiated in we 



address the problem of the 3PN dynamics of point-particles in harmonic coordinates. Earlier 
in [Q, our result has been already discussed and reported in the case of circular orbits. We 



find the presence of one (and only one) undetermined coefficient in the 3PN equations of 
motion, in agreement with [^^,^ . Recently, the physical equivalence between our result 
in harmonic coordinates and the result given by the ADM-Hamiltonian approach has been 
established p6| , |37[ . 

Another line of research, initiated by Chandrasekhar and collaborators [^,^], consists 
of working with continuous hydrodynamical fluids from the start, and derived the metric and 
equations of motion of an isolated fluid ball up to the 2.5PN order |]40|-^ (the derivation 
in the case of two fluid balls, in the limit of zero size of the bodies, being due to [p^ ,^). 
Our iteration of the gravitational fleld and equations of motion in the previous paper pTI is 



close to the latter line of work in the sense that it is based on the reduction of some general 
expressions of the post-Newtonian metric, initially valid for continuous fluids, to point-like 
particles. The choice of point-particles, adopted here as well, is motivated by the efficiency 
of the delta-functions in performing some complicated non-linear integrations. The price we 
have to pay is the necessity of a self-field regularization. We apply systematically in this 
paper the regularization of Hadamard, based on the concept of "partie finie" of singular func- 
tions and divergent integrals p4|-|46[]. This technique is indeed extensively used in this field 
29|j3T|,|32| . More precisely, we apply a variant of the Hadamard regularization, together 



with a theory of pseudo-functions and distributional derivatives, that is compatible with the 
Lorentzian structure of the gravitational field. All the details about this regularization can 
be found in We use notably a specific form of distributional stress-energy tensor 

based on "delta-pseudo- functions" (with support limited to the world lines of the particles). 
In a sense, these delta-pseudo-functions constitute some mathematically well-defined ver- 
sions of the so-called "good delta functions" introduced long ago by Infeld (see also an 
appendix in the book of Infeld and Plebanski |p(]|). 

Thus, we are using a formal regularization method, based on a clear mathematical frame- 
work 11^,48], but that we cannot justify physically (why should the compact objects be 
described by such delta-pseudo-function singularities?). Definitely, our main justification is 
that this method permits the derivation of a result in a consistent and well-defined way (i.e. 
all the difficult non-linear integrals at the 3PN order are computed unambiguously). Fur- 
thermore, we shall check that some different regularization prescriptions yield equations of 
motion that are physically the same, in the sense that they differ from each other by merely 
a coordinate transformation. Moreover, a justification a posteriori is that the end result 
owns all the physical properties that we expect the true equations of motion of compact 
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objects to obey. In particular, there is agreement with the known results at the previous 
post-Newtonian orders, we get the correct geodesic limit for the motion of a test particle 
in a Schwarzschild background, find that the 3PN equations of motion stay invariant under 
global Lorentz transformations, and obtain a conserved energy at 3PN (neglecting the ra- 
diation reaction). The investigation of the Lagrangian formulation of the equations is dealt 
with in a separate work fS^ . 

Ideally, one should perform, instead of a computation valid only for point particles (and 
necessitating a regularization) , a complete calculation in the case of extended bodies, i.e. 
taking into account the details of the internal structure of the bodies. By considering the 
limit where the radius of the two objects tend to zero, one should recover the same result as 
obtained by means of the point-mass regularization. This would demonstrate the suitability 
of the regularization. In fact, this program has been achieved at the 2PN order by Grishchuk 
and Kopejkin P7|J28| , who proved that the compactness parameters associated with each 
object disappear from the equations of motion, and obtained the same equations as in the 
case of point particles. At the 3PN order there is no such a proof that the method with 
extended bodies would give the same result as with point particles. 

The main problem is that from the 3PN level one cannot compute the most difficult of 
the non-linear integrals in closed form for two extended fluid bodies of finite radius (though 
these integrals could perhaps be obtained as power series valid when the two radius tend to 
zero). Presently the only approach which is able to overcome this problem is the one followed 
in this paper: namely, to model the source by delta functions and to use a regularization. 
The price we have to pay is the appearance of one physical undetermined coefficient at the 
3PN order. As a consequence, this method should be completed (hopefully in a future work) 
by the study of the limit relation of the point-particle result with the physical result valid 
for extended bodies in the limit of zero size. This study should probably give the value of 
the undetermined parameter left out by the regularization. 

The plan of this paper is the following. In Section |I|, we review some necessary tools 
concerning the regularization and the definition of the point-particle model. In Section |T|, 
we perform the post-Newtonian iteration of the field equations and write the 3PN metric in 
terms of some convenient non-linear potentials. Section |V| is devoted to the computation of 
the compact-support and quadratically non-linear parts of the potentials. The most difficult 
potentials, involving notably some non-compact cubic non-linearities at IPN, are obtained 
in Section 0. The so-called Leibniz and non-distributivity contributions to the equations 



of motion are derived in Section Finally, we present in Section [VI 1| the result for the 
compact binary's 3PN acceleration (in the case of general orbits) and the associated 3PN 
energy. 



II. HADAMARD REGULARIZATION 



In this section we present a short account about the regularization of Hadamard |44,45 



the associated generalized or pseudo-functions, and the choice of stress-energy tensor for 
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point-particles. We follow (and refer to) the detailed investigations in p7| , ^ . Consider the 
class T of functions -F(x) which are smooth (C°°) on deprived from two singular points 
yi and y2, around which they admit a power-like singular expansion of the type 



VnGN, F(x)= r?/„(ni) + o(r;^), (2.1) 

aQ<a<n 

and similarly for the other point 2. Here t\ = |x — yi| — 0, and the coefficients if a of 
the various powers of ri depend on the unit direction ni = (x — yi)/ri of approach to the 
singular point. The powers a of ri are real, range in discrete steps (i.e. a G (ai)^^^) and are 
bounded from below (oq < a). The coefficients ifa (and 2/a) for which a < are referred 
to as the singular coefficients of F. If F and G belong to so does the ordinary pointwise 
product FG, as well as the ordinary gradient diF. We define the Hadamard "partie finie" 
of F at the location of the singular point 1 as 



iF\ = J ^/o(ni) , (2.2) 

where dQi = dQ{ni) denotes the solid angle element centered on yi and of direction ni. 
Furthermore, the Hadamard partie finie (Pf) of the integral J d'^x F, which is in general 
divergent at the two singular points yi and y2 (we assume no divergence at infinity), is 
defined by 



Ffs,,s2 j (T^F = lim <i / d^^F 



3, 



+ ^T,y +47rlnf-Vr?F) +1^2!. (2.3) 

The first term integrates over a domain 'D{s) defined as to which the two spherical balls 
Ti < s and r2 < s of radius s and centered on the singularities are removed. The other terms, 
in which the value of a function at 1 takes the meaning ( p.2|) , are such that they cancel out 
the divergent part of the first term in the limit where s ^ (the symbol 1 2 means the 
same terms but corresponding to the other point 2). Note that the Hadamard partie finie 
depends on two strictly positive constants si and S2, associated with the logarithms in ( ^.31 ). 
See and Section |V| below for alternative expressions of the partie-finie integral. 

To any F G we associate a partie finie pseudo-function PfF defined as the linear form 
on JF given by the duality bracket. 



VG G ^ , < PfF, G >= Pf / d^x FG . 



(2.4) 



The pseudo-function PfF, when restricted to the set of smooth functions with compact 



support, is a distribution in the sense of Schwartz ||4^. The product of pseudo-functions 
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coincides with the ordinary pointwise product, namely PfF.PfG = Pi{FG). A particularly 
interesting pseudo-function, constructed in on the basis of the Riesz delta function 
is the delta-pseudo-function Pf^i, which plays the same role as the Dirac measure in 



distribution theory, in the sense that 



WF eJ^ , < PfSi, F >= Pf j rf^x 6iF = (F)^ 



(2.5) 



where [F)-^ is the partie finie of F as defined by ( |2.2|) . From the product of PfSi with any 
PfF we obtain the new pseudo-function Pi{F6i) which is such that 



VGg^, <Pf{F6i),G>= (FG)^ 



(2.6) 



Next, the spatial derivative of a pseudo-function of the type PfF, namely di{PfF), is 
treated as follows. Essentially, we require in ||4^ the so-called rule of integration by parts, 
namely that we are allowed to freely operate by parts any duality bracket, with the all- 
integrated ("surface") terms always zero exactly like in the case of non-singular functions. 
This requirement is motivated by our will that a computation involving singular functions 
be as much as possible the same as a computation valid for regular functions. Thus, 



VF,GgJ^, <d^{PfF),G >= - <di{PiG),F > 



(2.7) 



Furthermore, we assume that when all the singular coefficients of F vanish, the derivative of 
PfF reduces to the ordinary derivative, i.e. di{PiF) = P{{diF)- As a particular case, we see 
from these assumptions that the integral of a gradient is always zero: < di{PiF), 1 >= 0. 
Certainly this should be the case if we want to apply to the case of singular sources a formula 
which is defined modulo a total divergence for continuous sources. We have also at our 
disposal a distributional time derivative and the associated partial derivatives with respect 
to the points 1 and 2 (see Section IX in |^). The difference between the distributional 
derivative and the ordinary one gives the distributional terms D j [F] present in the derivative 
of F, 



a,(PfF) = Pf(9,F) + D,,[F]. (2.8) 
A simple solution of our basic relation (|2.7| ), denoted D^'^'^*[F] standing for the "particular" 



solution, was obtained in |47| as the following functional of the singular coefficients of F, 



1 x 1 X 



fc>0 



r7 1 



-2-fc 



(2.9) 
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where we assume for simplicity that the powers a in the expansion of F are relative integers, 
a G Z. (The sum over k is always finite.) The distributional term ( |2.9| ) is of the form Pf 
(plus 1 2). However, the particular solution ( |2.9| ) does not represent the most satisfying 
derivative operator acting on pseudo-functions. It is shown in that one can require also 
the rule of commutation of successive derivatives, which is not satisfied in general by ( ^.9p . 
Still we are motivated when asking for the commutation of derivatives that the properties of 
our distributional derivative be the closest possible to those of the ordinary derivative. The 
most general derivative operator satisfying the same properties as ( p.9|) and, in addition, the 
commutation of derivatives (Schwarz lemma) is given by 



+00 



D,[F] = 47r^Pf Q 



1=0 



1 1 . 



^i^i + Etx/V.^i +1-2 



(2.10) 



fc>0 



1 1 



where we denote by 1/^ the STF-harmonics of the expansion coefficient i/^, which is such 
that if a = ^i>Qf^iift (see |^ for details). A particularity of this derivative is that it 
depends on an arbitrary constant K through the /-dependent coefficient 



Ci = {l + 1) 



f2.111 



Both the derivative operators ( p.9|) and (|2.1CI| )- (p.ll|) represent some generalizations of the 
Schwartz distributional derivative that are appropriate to the singular functions of the 
class T. In was shown in Section VIII in that the distributional terms associated with 
the /th distributional derivative, i.e. Di[F] = S^PfF — PfS^F, where L = iii2 . . .ii denotes 
a multi-index composed of / indices, is given by 



^l[F] =J29n...,-Md^,^....HF] . (2.12) 



k=l 



Though this is not manifest on this formula, D2,[-F] in the case of the "correct" derivative 
( ^.1U| )-( P.11| ) is fully symmetric in the I indices forming L. Note that neither of the deriva- 



tives ( p.9|) and (|2.10|) satisfy the Leibniz rule for the derivation of a product. Rather, the 



investigation in |^ has suggested that, in order to construct a consistent theory (using 
the "ordinary" product for pseudo-functions), the Leibniz rule should in a sense be weak- 
ened, and replaced by the rule of integration by part (|2.7|), which is in fact nothing but 
an "integrated" version of the Leibniz rule. In this paper, we shall be careful about taking 
into account the violation of the Leibniz rule by the distributional derivative. We shall also 
investigate the fate of the constant K appearing in (p.ll|) when deriving the 3PN equations 
of motion. 

The Hadamard regularization [F)^ is defined by ( |2.2D in a preferred spatial hypersurface 
t =const of a coordinate system, and consequently is not a priori compatible with the global 
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Lorentz invariance of special relativity. If we restrict the coordinates to satisfy the usual 
harmonic gauge conditions, we introduce a preferred Minkowski metric, and thus we can 
view the gravitational field as a relativistic Lorentz tensor field in special relativity, that 
we certainly want to regularize in a Lorentz-invariant way. To achieve this we defined in 
^ a new regularization, denoted [F] ^ , by performing the Hadamard regularization within 



the spatial hypersurface which is geometrically orthogonal (in the Minkowskian sense) to 
the four- velocity of the particle. In a sense, the regularization [F]^ permits us to get rid 
of the anisotropic Lorentz contraction due to their motion when defining the point-masses. 
The Lorentzian regularization [F]-^ differs from the old one (F)^ by relativistic corrections 
of order 1/c^ at least. All the formulas for its computation are given in in the form 



of some infinite expansion series in the relativistic parameter The regularization [F]-^ 

plays a crucial role in the present computation, as it will be seen that the breakdown of 
the Lorentz invariance due to the old regularization (F)^ occurs precisely at the 3PN order 
in the equations of motion. Associated with the new regularization in we can define. 



exactly like in ( |2.5|) , a "Lorentzian" delta-pseudo-function PfAi which when applied on any 
F gives [F]-^^. More generally we have, similarly to (|2.6| ), 

VGgJ^, < Pf(FAi),G>= [FG]^ . (2.13) 

Notice that as a general rule we are not allowed to replace F in the pseudo-function Pf (FAi) 
by its regularized value, i.e. Pf(FAi) 7^ [F]-^ PfAi. This is a consequence of the "non- 
distributivity" of the Hadamard partie finie with respect to the multiplication, i.e. [FG]j^ 7^ 
[F]-^ [G]-^. In this paper, we shall (heuristically) model the compact objects by point-particles, 
and in order to describe those point-particles we shall use a particular representation of the 
stress-energy tensor which has been derived in Section V of [EH] on the basis of an action 



principle compatible with the Lorentzian regularization [F]-^. The proposal made in |^ is 
that 

T/^- = m,c^jt^=Pf ( 4^] +1^2. (2.14) 

Most importantly about this expression are the facts that (i) [gpa]i within the first factor 
means the Lorentzian regularization of the metric in the previous sense, (ii) the pseudo- 
function Pf ^-^Aij is of the type Pf(FAi) which is defined by (|2.13|) . We denote by mi 
the (constant) mass of the particle 1, by yi(t) its trajectory parametrized by the harmonic- 
coordinate time t, and by vi(t) = dyi/dt the coordinate velocity [with Vi = (c, vi)]. In the 
next section, we look for solutions in the form of post-Newtonian expansions of the Einstein 
field equations having the latter stress-energy tensor as a matter source. 

III. THE THIRD POST-NEWTONIAN METRIC 
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A. The Einstein field equations 



We base our investigation on a system of harmonic coordinates x° = ct, (a;*) = x, 
since such coordinates are especially well-suited to a post-Newtonian (or post-Minkowskian) 
iteration of the field equations. We define the gravitational perturbation h'^'^ associated with 
the "gothic metric" as 



h'''' = V^g'''-v''\ (3.1) 

with g'^'^ and g being the inverse and the determinant of the covariant metric gfj,^, and 
where t]^'^ = diag(— 1, 1, 1, 1) denotes an auxiliary Minkowski metric. Under the condition 
of harmonic coordinates, 



d^h^'' = , (3.2) 
the Einstein field equations take the form 

nW"" = ^^|^|T^^ + A''^ , (3.3) 

where □ = rj^^d^d^, denotes the flat d'Alembertian operator, where T^'^ is the matter stress- 
energy tensor defined in our case of point-particle binaries by ( |2.14| ), and where A^'^ is the 
gravitational source term. Using the integral of the retarded potentials given by 

n-i^. ,^ / d'^' r(x^t-|x-xVc) 

a r(x,t) = / — ■ , (3.4) 

J — 47r |x — x'l 

we can also re- write the solution of the field equations ( |3.3|) , under a condition of no-incoming 
radiation, under the form 



C4 



(3.5) 



The gravitational source term A'^'^ is related to the Landau-Lifchitz pseudo-tensor by 

1 Rjrf^ 

A^'' = ^\gKl + d,h^^d,h^P - h^-d.^h'^-' , (3.6) 

and can be expanded as an infinite non-linear series in h and its first and second space-time 
derivatives; in this paper we need only the non-linear terms up to the quartic {h^ or G^) 
level, viz 
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p^i^v ^ Ni'^^h, h) + M^'^(/i, h, h) + L>'''{h, h, h, h) + OQv' 



(3.7) 



where the quadratic non-hnearity N^^ , the cubic one M^^ and the quartic L^^ are exphcitly 
given by 



1 



1 
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pit-crr'- 



hP'^f^dPhprd'hl + drh^;d^hl - dph'^d^h"^^ 

+ h"" - -drhp^d'^h'^ + \dphdPh + \dpKrd^h'^ 
4 8 2 



d'^Krdph"^ - 2d^Kpdph''^ - -d'^hdph 



1 



1 

.p..^a-i - ^ 
1 



1 

4' 



- ^h^^dphra&'hl + -h^%KdPhxa 

L^- = - Ik'^^hpadrhP^dxh"^ - ^h'^^hpadPhrxd'^h''^ + lh'"'hpad''hd^h 
2 4 8 



1 

2' 



+ ^hi^'^hp^drhld^h"^ - ^h^^^hp^drhP^d^h + hpxKdrh^Pd^h' 



2hpxh%h"^^d''^h''^ + Upxh^^d^^hP^d^h"^ - ^hpxh^ad^^h'"'d''^h 



- hpahrxd'^h^^'d^h'"' + ^hpahrxd'^h'^&'h''^ - ^hpahrxd''h'"'d''h^^ 

+ 2KrhP^^dxh''^^dp}f^ - ^KM^&'^ldph^^ + ^-harh^^^d^'^''^ dph 

+ "^-KrhF^^&'^hdph"^ + -h^'Ph'^'dphrxdah^^ - ]-hi'PK"'dphd^h 
2 2 ^ 



+ T] 



^hp^hldrhP^dxh'^- - hp^hldxKd^h''^ + ^hp^hldrhP'^d^h 



- Ihp^Kxd^hP^d'^h''^ + hp^hrxd^h'^'d-'h^^ + hp^hrxd^Kd'^h 

4 o Z 

- ^hp^hrxdPh^^&'h 



Att 



(3.8a) 



(3.8b) 



(3.8c) 



All indices are lowered and raised with the Minkowski metric rj^iy] h — rj^'^hniy] the parenthesis 
around indices indicate the symmctrization. 

To describe the matter source we find convenient to introduce the density of mass a, of 
current cTj and of stress aij defined by 



(3.9a) 
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(T,;C 



Oi 



r 

rpij 



(3.9b) 
(3.9c) 



— SijT^^). These definitions are such that a, ai and aij admit a finite non-zero 



(where T 

hmit when c +00 (since T^^ has the dimension of an energy density). In the case of our 
model of point-particles [stress-energy tensor given by (p.l4|) ], we obtain 



a(x,t) = Pf(/iiAi) + 1^2, 
ai(x,t) = Pf(/iit;iAi) + l ^2, 



(3.10a) 
(3.10b) 
(3.10c) 



where Ai = A[x — yi(t)], and where /xi and /ii represent some effective masses defined by 

ITLiC 1 



/xi(x,t) 



/ii(x,t) = /ii(x,t) 



(3.11a) 
(3.11b) 



Note that /xi and jli depend on time and space. Indeed, while the first factor in (p.llaj) 
is clearly a mere function of time through the values of the positions and velocities of the 
particles at the instant t, the second factor (—(7)^^^^ is evaluated at the /le/d point t, x instead 
of the source point t, yi. From the non-distributivity of the Hadamard regularization, one 
is not allowed to replace {—g)~^/'^ by its (regularized) value at the point 1, even though it 
is multiplied by a delta-pseudo-function at 1. 



B. The 3PN iteration of the metric 

In what follows we sketch the main steps of our iteration of the Einstein field equations 
( p.3|) -( p78D generated by two particles at the 3PN order. For more clarity in the presentation, 
we reason by induction over the post-Newtonian order n. However, we do not have proved 
the validity of this method to any order n; simply we applied the method outlined below to 
construct the metric at the 3PN order. 

(I) Suppose by induction over n that we have succeeded in obtaining some approximate 
post-Newtonian metric coefficients h'^2n-2] ' ^^^^ previous coefficients /ij^j for any m 

such that 2 < m < 2?2 — 2, which approach the true metric modulo a small post-Newtonian 
remainder. 



h^" = h^L_,^+0{2n-l), (3.12) 

with the notation 0{2n — 1) = (9(l/c^"~^), here and elsewhere, for the post-Newtonian 
error terms. We assume that the /^Ci's are at once some explicit functions of the field point 
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X and functionals of the two trajectories yi(t), y2(t) and velocities vi(t), V2(t). Since the 
matter source of the field equations is made of delta-pseudo-functions, the metric coefficients 
become singular at the location of the particles (indeed, this is already true at the Newtonian 
order). As a matter of fact, we assume for the present iteration that. 



Vm < 2n - 2 



Im] 



(3.13) 



where is the class of functions considered in ||47| , |48[| and Section 0. This is not a completely 
rigorous assumption because of the presence of logarithms in the expansions around the 
singularities; but we shall see that this assumption is justified at the 3PN order where one 
can consider these logarithms as mere constants. 

(II) Consider for simplicity the combination |(/i°° + /i**) only, for which we need the maximal 
post-Newtonian precision since it is directly connected to Qqq. The structure of the Einstein 
field equations (|3.3| ), containing notably the gravitational source term (|3.7|) - (|378|) , reads as 



□ 



SttG 



g\a + h . . . hdhdh , 



(3.14) 



where a is given by ( p.9a{ ), where the sum runs over non-linearities and the two partial 
derivatives dd have to be distributed among the /I's (with double derivatives allowed in the 
quadratic term). In order to obtain an equation valid at the next post-Newtonian order n, 
we replace the approximate metric (|3.12| ) into the right-hand side of ( p.l4|) . Furthermore, 
we replace the partial derivatives dd in (|3.14|) by the distributional derivatives ( p.8|) [we 
shall discuss the effect of using either the particular derivative ( |2.9D of the more correct one 
( p.lOj )]. Using also the density of particles in the form ( p.lOa|) , we get 



□ 



|^(pf|(?|^lAi+Pf|(?|/i2A2 



+ 

mi,...,mp<2n—2 

+ 0{2n + l) 



/i[mi] • • • h[mp_^]d (Pf 8 (Pf 



[2n] 



(3.15) 



where the h^mi], ■ ■ ■ , ^[mp] (with 2 < p < n) denote the metric coefficients known from the 
previous iterations, and where as indicated by the label [2n] a truncation up to the post- 
Newtonian order l/c^^ is understood. At this stage, any subsequent transformation of the 
right-hand side must be done using the rules for handling the pseudo-functions and their 



derivatives [^ . 

(III) We integrate the latter equation by means of the retarded integral given by 



Uii C RttC / 

— = I [pf \g\f^,A^ + Pf \g\fi2A2 



12 



+ 



h[m,] ■ ■ ■ h[m,_,]d (Pf d (Pf hlm^]) 



[2n] 



+ 0{2n+l) . 



(3.16) 



This defines the solution to the nth order, and so, by recursion, to any order (in principle). 
The partie-finie symbols Pf take care of the divergences of the retarded integral at the 
locations of the particles; that is, the retarded integral is considered as a partie-finie integral 
in the sense of [^]. More precisely, the retardations in ( p.l6| ) are expanded to the nPN 
order and the resulting Poisson-like integrals computed using the duality brackets in the 
way specified by Section V in [^. Actually, the Poisson-like integrals, which have a non- 
compact support, become rapidly divergent at infinity when n increases, and the correct 
solution we use is not the Poisson-like integral but is obtained by a matching of the inner 
metric to the multipole expansion of the exterior field. So, in fact. 



□^ipfF (x',t) 



2n 



Air ^ k\& 

k=0 



< 



PfF(x,t) 



^ match 



+ C(2n + 1) 



(3.17) 



where the subscript "match" refers to the matching process that is described in Section IV 
in the case of the 3PN order. Notice that the time-derivatives {d/dt)^ resulting from the 
Taylor expansion of the retardations are distributional derivatives and therefore can be put 
outside the duality bracket (see Section IX in [0). Thus, equivalently. 



(□^^PfF)(x',t) = - 



2" 



k=0 



< PfF(x,t), |x 



^ match 



+ 0{2n + 1) . 



(3.18) 



(IV) Once the solution ( p.l6| ) to the nth post-Newtonian order is in hands we perform 
many simplifications of the expression, following the rules of application of the distributional 
derivative. In particular, we find very useful to use the fact that a double gradient can be 
re-expressed in terms of d'Alembertians as 



d^Fd^'G = - [□ (FG) - FUG - GUF] 



(3.19) 



which implies that the retarded integral reads as 



= ]^FG - in^^i [FUG + GUF] . (3.20) 

The first term is "all-integrated", while the second term, in which one can replace the 
d'Alembertians by their corresponding sources, brings in general many interesting cancella- 
tions with other terms. Unfortunately, the formula (|3.19| ) is valid only in an ordinary sense 
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but not in the distributional sense, because the distributional derivative does not satisfy, in 
general, the Leibniz rule. Thus, in general. 



a^(PfF)a''(PfG) ^ - [□ (PfFG) - Fn(PfG') - G'n(PfF)] . (3.21) 

Nevertheless, the strategy we have chosen to follow in this paper is to take advantage of the 
many simplifications brought about by the latter process, at the price of introducing some 
extra terms (named "Leibniz") accounting for the violation of the Leibniz rule. This means 
that we shall write, similarly to (|3.20|) , 

|a,(PfF)a^pfG)| = iptre - iFD(pfG) + GD(PfF)i + .^„,.r . (3.22) 

where the Leibniz term is given by 



5 



Leibniz 



T = □ 



7^ 



9^(PfF)9^(PfG') 



in(PfFG) + ^Fn(PfG) + ^Gn(PfF) 



(3.23) 



Obviously the Leibniz term depends only on the purely distributional part of the derivative. 
See the Appendix ^ for the complete list of the Leibniz terms. As it will turn out these 
terms are not too difficult to compute, and, of course, arise precisely at the 3PN order. They 
give a contribution to the metric and the equations of motion that we shall be able to check 
from the requirement of Lorentz invariance (see Section IVT]) . 



C. The 3PN non-linear potentials 

The post-Newtonian iteration sketched in the previous subsection is implemented to the 
3PN order. The computation is long but straightforward. After the simplification process 
described above we find that the metric is parametrized by certain non-linear potentials, 
which do not carry a physical signification by themselves, but turn out to be useful in the 
present computation. The 3PN metric reads as 



9oo 



gij 



-1 + ^v^ 



+- 



32 
4 



T 



-VX + RiV, - -VViVi - —V^ 
2 * * 2 48 



+ 0(10) 



rRi 



16 



Y^ + \w.,V, + ^ ) + 0(9) 



^1] 



2 2.8 

& & & 



X + + 



+- 



16 



VWij - V,Vj + C(8 



(3.24a) 
(3.24b) 



(3.24c) 
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We recall our notation for the small post-Newtonian remainders: 0{n) = 0{l/c^). The 
various post-Newtonian orders are parametrized by some potentials which are defined by 
means of the retarded integral ( |3.4| ). At the "Newtonian" and IPN orders we pose 

(3.25a) 
(3.25b) 



V = □^i[-47rG(T] 
Vi = n:^'[-AnGa,] 



in which the source densities were defined by ( p.9|) . Next, at the 2PN order, we define 



X = l-iTiGVaa + Wi,d,,V + 2VAdiV + Vd^V + -{dtVf - IdiVp^V^ , (3.26a) 



R, = □^i[-47rG(y(Ti - Via) - 2dkVdiVk - -dtVdiV] 



(3.26b) 
(3.26c) 



Finally, at the 3PN order, we have 



T — D"^ 
^ — '-'7^ 



-diV.dtW,, + VVA^^V + 2V,d,Vd,V + -VAVdy + -V'dtV 



+ 5: 



Leibniz J 



(3.27a) 



Y- — 

ii — 



-AtiG ( -aki - aVVi + ]^akWik + ^cr^kYk + ^a^^V^i ) + WkidkiV^ - dtWtkdkV 



+d,WMdkVi - dkWadiVk - 2dkVdA - - ^-VdtVdy 



-2VdkVdkV, + Vd^V, + 2VkdkdtV, 



~l~ '^Leibniz ^ 



(3.27b) 



-4nGV{a,, - 6,,akk) - 2d^^^dtV,, + d^Vkd^Vk + dkV.dkV^ - 2d>^AV, 



i) 



-SijdkVm{dkVm - d^Vk) - -Sij{dtVf 



+ 6: 



Leibniz ^ij 



(3.27c) 



Note the presence in the 3PN potentials of the Leibniz contributions described in the previous 
subsection, which are due to the simplifications we did to arrive at these relatively simple 
expressions (with respect to what could be expected at the high 3PN order). The Leibniz 
contributions will be computed in Section 0. Of course, in the case where the matter source 
is continuous — an hydrodynamical fluid for instance — , the 3PN metric ( |3.24| ) and all the 
expressions of non-linear potentials are valid with simply the Leibniz contributions set to 
zero. 

The potentials (p.25|) - (|3.27|) are connected by the following approximate post-Newtonian 
differential identities (equivalent to the condition of harmonic coordinates at the 3PN order): 



dt<V + 



Wkk + 2V' 



+ 



X + ^Zkk + \vWuk + 



15 



+9, <^ V, + 



R^ + VV, 



+ 



- \^^^,V, + + VRi + r^V, 







1 


Zij - 







= 0(6) , (3.28a) 
C(4) . (3.28b) 



We shall check that the (regularized) potentials we compute satisfy these identities. They 
are in turn respectively equivalent to the equation of continuity at the 2PN order and the 
equation of motion at the IPN order: 



a\ 1 + 



dt 



CJi 1 



a,\ 1 



21^,;, 



i {dta,, - adtV) - i {aV,d,V + a.^d^Vk) + 0(6) 



era 1 + 



AV 



adiV + ^ [adtVi + aj {djV, - diVj)] + 0(4) 



(3.29a) 



(3.29b) 



D. Computing the equations of motion 

The equations of motion of the particle 1 are deduced from the covariant conservation 
of the stress-energy tensor of the particles, 



V^T"'' = 



(3.30) 



where T'^'^ is given by the definite expression ( p.l4| ) made of the delta-pseudo-functions 
defined in |^,^. It is shown in Section V of that by integrating ( p.30| ) over a volume 
surrounding the particle 1 (and only 1), i.e. by constructing the duality bracket of ( p.3U[ ) 
with the characteristic function of that volume, we obtain the equations of motion of the 
particle 1 in the form 




"1 "1 

per] 1 c2 



1 [d\gf^u]iViVi 

[y peril c2 



(3.31) 



These equations of motion take the same form as the geodesic equations for a test particle 
moving on a smooth background, but with the role of the background metric played by 
the true metric generated by the two bodies and regularized according to the Lorentzian 
prescription [^]. 

In this paper we compute the spatial acceleration of body 1, which corresponds to the 
equation with spatial index A = i in (|3.31|) ; we do not consider the energy which would 
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be given by the equation with time index A = 0. Indeed, the energy of the binary system 
will be determined directly from the (fully order-reduced) acceleration. From ( |3.31| ) we can 
write the equations into the form 



f^F;, (3.32) 

where the "linear momentum density" and "force density" are given by 

Fl = i J^£dl±L . (3.33b) 



2 



The expressions of both PI and Fl in terms of the non-linear potentials follow from insertion 
of the 3PN metric coefficients (|3.24|) . We obtain some complicated sums of products of 



potentials which are regularized at the point 1 following the prescription [F]-^. Since the 
computation will turn out to be quite involved, we decide to adopt the following "step-by- 
step" strategy: 

(A) We compute, in Sections |I^ and 0, all the needed individual potentials and their gradi- 



ents at the point 1 following the non-Lorentzian regularization (-F)^; for instance we obtain 
(diV)^ at the 3PN order, (V)-^ at the 2PN order, (diX)-^ at the IPN order, (dif)^ at the 
Newtonian order, and so on. (Because of the length of the formulas, and since the results 
for each of these individual regularized potentials are only intermediate, we shall not give 
them in this paper; see the appendices of [0] for complete expressions.) 

(B) We add up the corrections brought about by the Lorentzian regularization [F]^ with 
respect to (F)^. We find, at the end of Section 0, that the only effect of the new regu- 
larization at the 3PN order, when computing the values of potentials at 1 (but the new 
regularization affects also the corrections due to the non-distributivity), is a crucial IPN 
correction arising from the so-called "cubic non-compact" part of X; that is, we find 

(C) We replace all the individually regularized potentials [F]-^ and their gradients into 
the equations of motion (p.32|) - (|3.33|) which would be obtained while supposing that the 
Hadamard regularization is ^^distributive" with respect to the multiplication, i.e. suppos- 
ing incorrectly that we are allowed to write everywhere [-FG]^ = [-Fj-j^lG]^. In doing this 
we obtain what we call the "distributive" parts of the linear momentum and force densi- 
ties ( p.33|) , namely (-Pi)distr (-^i)distr- (Other types of non-distributivity arising in the 
potentials themselves are discussed in Section |V[) 

(D) Finally, we compute separately, in Section VI, the corrections due to the non- 
distributivity, i.e. the differences -Pi — (-Pi) distr -^i ~ (-^i) distr- Note that these corrections 
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reflect quantitatively the specific form that we have adopted for the stress-energy tensor of 
point-particles ( |2.14| ). Had we used another stress-energy tensor, for instance by replacing 
incorrectly Pf(;^^Ai) by [;^^]-,^PfAi inside ( |2.14| ), we would have obtained a different non- 
distributivity, and thereby some different equations of motion. Note also that thanks to the 
new regularization [F]-^ the corrections due to the non-distributivity do not alter the Lorentz 
invariance of the equations of motion. At last, we find the 3PN acceleration of body 1 as 



a\ = Fl - ^^[Pl - v\) . (3.34) 

We report now the expressions of the distributive parts of the linear momentum and 
force densities as straightforwardly obtained by substitution of the 3PN metric (|3.24| ). The 
expressions of the correcting terms due to the non-distributivity (i.e. [FG]^ ^ [-^]i[^]i) ^.re 
relegated to Section Vl, where it is seen that they contribute only at the 3PN order. 

(A^)d.st. = + ^ (^^^i + m.v\ - 



2' 

9, 



- m,\A< + \^fA + m,\mA 

-m^3\\nx - 8^1 - 16[F.]0 + 0(8) , (3.35a) 



+ \, i^\mA - mv,]AA + livumA 



+2[d..w,u]AA - mAmAi - mi[dv,]A 

-mR3]Ai + \[v]l[d,v], + s[v,Ud,v,],+A[d,x], 

1/11 3 . 4.Q 

-nvumAAi - im,m]AAi - mR3]AAi 

27 

+-[v]l[dy],vi +i2[v,],[d.v,],vl + mMmAA 
+Q[v\[d,w,u\vA + mx]A + mZjk]AA - '^^\y3W]A^v]Al 
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-A[X\[diV\ - 4[y]J9,X], + imT\) + 0(8) . (3.35b) 



Recall that it is supposed that all the accelerations appearing in the potentials have been 
order-reduced by means of the equations of motion. Notably, during the reduction of the 
"Newtonian" term [^jV^]^ in Section we shall need the equations of motion to the 2PN 
order. Furthermore, we see from ( p.35a| ) that when computing the time-derivative of Pi 
we meet an acceleration at IPN which is thus also to be replaced by the 2PN equations 
of motion. We recall here that the latter 2PN (or, rather, 2.5PN) equations in harmonic 



coordinates are 24-26,31 



dv\ 



Gm-2 



vl^ [4:{ni2Vi) - 3(ni2t;2)] 



12 2 ^2 

' 12 ' 12'- 

+n*i2 -^1 - 2t^2 + 4:{viV2) + ^{ni2V2f + 5^^^ + 4- 

2 ri2 ri2 



. Grrii Gm2 



+ 



3 9 
-2v2 + 4:vl{viV2) - 2{viV2f + -vfini2V2y + -jVl{nx2V2f 



15 



-'0{VxV2){nx2V2) - —{nx2V2) 



+ 



Cm I 
ri2 



-^^1 + ^^2 - 2(^1^2) 



39 17 

+ Y(?^i2^^i)^ - 39(ni2t;i)(ni2t;2) + y(?^12^^2)^ 

+ [4t;2 _ 8{viV2) + 2{ni2ViY - 4:{ni2Vi){ni2V2) - 6(ni2t;2)^] 
12 



+ - 



12 



■ 57 2 Q 2 69 
——rrii — 9m2 — —mim2 



+ 



Gm2 
'12'- 



Vi2< vl{ni2V2) + '^vl{ni2Vi) - ^vl{ni2V2) - '^{viV2){ni2Vi) 



9 



+4(fif2)(ni2f2) - 6(ni2t;i)(ni2^^2)^ + 2(^12^2)^ 



+ 



Grrii 
ri2 
AG'^mim2 



63, ,55, , 

-^(«12^^l) + ^(^12t^2) 



+ [-2(ni2t;i) - 2(ni2t^2)] 

^"12 



'^l2('^12fl2) 



Grrii 52 Gm2 n 
-6 + + 3v^ 



ri2 



3 ri2 



12 
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Grrii Gm2 



12 



— V 



ri2 



ri2 



12 



■+0(6) 



(3.36) 



Unavoidably, because of the proliferation of possible terms, the equations of motion at the 
next 3PN order are even much longer [see ( [7.16| ) below]. 



IV. COMPACT SUPPORT AND QUADRATIC POTENTIALS 

All the potentials that enter the linear momentum (-Pi)distr force density (-^i)(jistr 

are computed at the point 1 by means of the Lorentzian regularization [F]-^. However, we 
shall first determine their Hadamard partie finie in the usual sense (-F)^, i.e. by approaching 
the singularity in the spatial slice t = const. The difference between the two regularization 
processes does not affect any compact or quadratic potentials. 



A. Iterative computation of compact support potentials 



In this paragraph, we are interested in the compact terms involved in the equation of 
motion [see ( |3.35aD and ( |3.35b| )]. According to our previous remark, it is sufficient to evaluate 
them with the classical Hadamard prescription. We need (diV)^ up to 3PN order; (V^)^, 
{V^)„ and {d,Vj), at 2PN; {wff^)„ {^^w'j^^),, (rP),, (mf^), and (9,X(c))^ at IPN. 



^(C)n 



The remaining contributions are Newtonian: (X(C))^, (^iT^^))^, (>"i Oi and ((9,F/ j-^. 



We 



follow the same classification and nomenclature concerning the various parts of potentials 
— compact, non-compact, etc. — as in Section II of The compact (C) potentials are 
generated by sources with (spatially compact) support limited to the particles; for instance, 
V^(C) = V, and, from {^W^, 



fic) 



□ 



n 



'1 1 
-47rG ( -(TijWij + ^V^CTii + aViVi 



(4.1) 



Thus, by definition, the source ^(x, t) of each compact potential P^^-* is made of Dirac 
pseudo-functions, multiplied by some functions of the class JF: 

5(x,t) = Pf(FAi) + Pf(G'A2) , 

with F, G ^ J-". As a result, it is in general possible to find an explicit expression of P(C) 
over the whole space (for any x). Besides, the expansion under the integration symbol of the 
retardation of 5'(x',t— |x — x'|/c) as c goes to infinity is perfectly licit, because the integrand 
has a compact support: 



p(C) 



A-K ^ n!c" 

n=0 



/in-l 



S(x', t) 
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1 +°° / \n 



47r — ' n\c 

?i=0 



X 



[G|x 



i\n-\ 



(4^2) 



The sources ^(x, t) are supposed to be known at the current order. This imphes to proceed 
iteratively as explained in [|31]]. The reader is referred to this paper for more details. In 
short, we start from the V and Vi potentials, whose sources do not depend on any other 
ones at the lowest order. Indeed, we have = — 47rG'Pf (/tiAi) + 1^2 (and similarly for 
Vi), where fli = mi + 0{2), as it follows from insertion of the "Newtonian" metric into the 
definition ( p.llD of the effective mass. Hence, 



V = G 



dh 



Pf [/iiAi]--atPf[|x-x'|/iiAi] 



0{2) 



Gmi 



—dtmi + 1 
c 



2 + C(2) 



Gmi Gm2 



ri 



r2 



+ 0{2) 



(4.3) 



To obtain the regularized metric (at the location of the first body, say), we need the partie 
finie of the potential V at point 1, (V)^. Since we use Hadamard regularization, it is 
simply given by the value of its non-singular part when x = yi. Here, we find (V)^ = 
Gm2/ri2 + 0{2), with the notation = |yi — y2|. 

The computation of more complicated compact terms necessitates the knowledge of the 
effective masses /ii and fli beyond the Newtonian approximation. By substituting to g^^, 
the explicit 3PN expression (|3.24|) for the metric in the equations (|3.11|) , we get the general 
forms of both effective masses. As an example, fii at 2PN order reads 



mi 



1 + 



+ 0{6) 



-2V+[V], + 



■2Wu + 2V^-2V[Vl + -[V]l 



[V\ - 3Vvl + l[V]y, - A[Vlv[ + ^-vt 



(4.4) 



where we are careful at distinguishing the potentials computed at the field point x from those 
computed at the source point yi, and where we take into account the non-distributivity of 
the regularization (/ii differs only by some numerical coefficients). Thus, as emphasized in 



Section |ITI| , /ii and fli are functions of time and space. Replacing them by the regularized 
quantities {fii)^, (/ii)i (and 1 2) is definitely forbidden because, on one side, the partie 
finie is not distributive, and, on the other side, the usual Hadamard regularization does not 
coincide with the Lorentzian one. However, this replacement does not modify any compact 
potentials, with the notable exception of the 3PN contributions in V (see below). It is thus 
convenient to pose: 

Kiistr = □7^' [-AnG{fii),Pf6i-AnG{fi2)2PiS2] , 

and to calculate Vdistr and V^— Vdistr separately. In the other compact sources, we shall employ 
{fj'i)i, {p'i)i, etc. instead of /ii and fii for practical calculations at the 3PN approximation. 
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Furthermore, in all the compact terms, the action of the Lorentzian delta-pseudo-functions 
Pf Ai and PfA2, remarkably, reduces to the one of Pf(5i and Pf(52. From what precedes, it 
becomes obvious that, after the evaluation of (/ii)^ or (fli)^ at a given post-Newtonian order 
n, we can determine all the potentials to the precision As all the terms involving 

the retarded potentials in /ii appear already with a factor at least, we are then in a 
position to compute the right-hand-side of the equation The process is initiated by 

the computation of the Newtonian value of V as presented above. Most of the quantities 
needed to get {fli)^ at the 3PN order are obtained in Finally, the regularized value of 
yUi at point 1 is 



mi 



1 + ^ 



Gm2 3 r, 

h -f 1 

ri2 2 ' 



+ -T 



ri2 
3 Gm2 



+ 



2 ri2 



1 



G^mim2 / 3 Grrii 39 Gm2 
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35 



41 



8G^mim2 , 
H — {{nuvi, 



IGmi 

- irii2V2)) 



12 



41 
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2 ri2 4 ri2 

+^(ni2t;i)(ni2f2) - ^(^12^2)^ 
4 o 



+ y^l - ^K^2) + y^^2 - g(^12^^l) 



+ 



+ 



Ghr4 

rl2 
Gm2 

ri2 



-vl - {ni2V2f - ^{viV2) - {ni2Vi){ni2V2) + -{ni2Vif - 
2 2 2 ri2 



3 Gm2 15 



+ 
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'2{viV2f + hvlvl - 10(i;it;2)^^i H vl + 2w 



-(ni2t;2) ^2 - 0(^121^2) 

Z o 



11 



0(7). 



(4.5) 



-4i;2(fit;2) + 2(ni2t;2)^(fi^^2) - -^^i(^i2^^2)^y . ^g- 

In our notation, two vectors Vi, V2 between brackets represent the scalar product: (^1^2) = 
v\v2^ vl = v{v{. We recall that it is important to keep the grouping of factors imposed by 
the regularization in products of potentials. For instance: (VWij)^ 7^ (^)i(^ii)i- 

Among the compact potentials, the 3PN value of V is certainly the most difficult one to 
obtain, since the other quantities require only lower orders in powers of 1/c. We shall focus 
on Vdistr to illustrate the method we have followed. The difference V — Vdist will be handled 
in the next subsection. We begin with specializing the general formula for V to the case of 

Vdistr : 

Gi^jll) ^ G ^ , ^ ^ G ^0 r / ~ \ n Gr%'^r/~\ 91 G 



distr 



ri 



24c4 



5f [(/ii)i^?] 



(4.6) 



Since the Schwarzschild mass mi is constant, dtjli/c is of order 0(3) and does not contribute 
at the IPN level. For convenience, we shall introduce some special notation for the terms 
that occur at this approximation; we pose: 



U 



ri 



r2 



and X = G(/ii)^ri + G'(/i2)2'^2 • 
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Actually the potentials are to be considered as pseudo-functions and it is understood that 
there is a symbol Pf in front of them. Notably, the time derivatives appearing in ( [4.6| ) are 
distributional. The regularized effective mass as well as the distance to the first body 
ri depend on time through the trajectories yi,2(^) and velocities vi 2(t). We explicit the 
time-differentiations and obtain, at IPN, 



K,i.t. = U+ ^d^X + 0(3) 



Gnii 



1 + 



01712 

ri2 



+ ^ {-a\d,n + vXdln) + 1 



(4.7) 
2 + 0(3) . 



The accelerations are order-reduced by means of the equations of motion at previous orders. 
Notably, for computing the IPN term ■^d'^Pi[{jli)^ri] at relative order 3PN, we need the 
2PN acceleration given by ( |3.36|) . Once we have got Vdistr all over the space, the last stage 
consists of regularizing it, as well as its gradient, at x = yi using the Hadamard partie 
finie. Now, Vdistr can be divided into two parts, Vdistr n and Vdistr r2 corresponding to the 
sources —4:7iG{fli)^Pf6i and — 47rG'(/i2)2Pf'52 respectively. The first part, Vdistr ri, depends 
on X through ri only, and contains many terms that are either singular or vanish when 
X — s> yi, giving no contribution to the partie finie; on the opposite, the smooth terms with 
odd 1/c-power factors in ( |4.6| ) generally contribute. The part Vdistr r2 does not necessitate 
any regularization since it is already regular in the neighbourhood of x = yi. 

The remaining potentials are determined in the same way. However, we have to apply 
properly the formalism developed in [0. In particular: 

(1) The regularized value of some potential P*^*"^ is the partie finie of P*^*"-* computed initially 
outside the singularity. In the case where P*-*"-* is the Poisson integral of a compact source 
Pf(P(5i), with P G J^, we must take care [0 that 



[We generally do not write the dependence of the integrand on the integration variable, as it 
is evident from the context; thus, F6i is computed at point x' in the intermediate expression 
of ( |4.8|) and at point x in the last one.] 

(2) If P is not regular at point 1, we generally have Pf(P5i) ^ (P)^Pf5i, even when both 
members act, in the sense of pseudo-functions, on smooth test functions. This distinction 
is crucial, for instance, in the determination of T^'"^ at Newtonian order. Indeed, one of its 
contributions [first term in ( [4.1| )], denoted by T^'^^\ reads as: 

r(ci) = Pf I + Oil) 

= lGm,v[vi{j:^-^] +1^2 + 0(1) 
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which is different from Gmivlvl{Wij)^/4ri + l ^ 2+0{l) . Had we used the latter expression 
instead of T^^^\ we would have obtained a different potential Tdistr! this would have been 
correct if the partie finie operation had been "distributive" (see Section |ID, but we have 
actually 



distr 



1 2 1 

vi H 1 

240 ^ 80 



nivi] 



+ 1^2 



Notice that the latter expression is not Galilean-invariant by itself, and therefore will be 
checked later when verifying that the final equations of motion stay invariant under Lorentz 
transformations. 



B. Non-distributivity in the potential V 



We call non-distributivity in the potential V that contribution which arises because the 
coefficient of the delta-pseudo-function Pf Ai in the matter stress-energy tensor ( p.l4| ) is a 
function not only on time but also on space through the factor 1/ ^—g. It will turn out that 
this contribution is purely of order 3PN. A related contribution, due to the non-distributivity 
in T, has just been computed in the previous subsection. The potential V is generated by 
the source density ct(x, t) = /iiA[x — yi(t)] + 1 2, where /2i is a function of space-time 
given explicitly by 



rriic ( 1 + - 

The first factor is a function of time, and the second one depends on both time and space 
(non-distributivity). The potential V is given by the retarded integral (|3.25a|) , whose retar- 
dations we expand up to any post-Newtonian order: 

r(x,t) = G5^i^9r / t/v |x-xr~v(x',t) . (4.10) 



n=0 



(Actually we shall see that an expansion to the IPN order is sufficient for our purpose.) 
We insert into that expression the source density a, use the definition of the delta-pseudo- 
function Pf(FAi) given by ( |2.13| ), and arrive at 



+00 



n=0 



/in-l 



(4.11) 



Here, the square brackets refer to the Lorentzian Hadamard regularization when x' — >■ yi. 
Using a multipolar expansion, we obtain immediately these brackets as 
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/ii(x',t)|x-x' 



/in-l 



1=0 



r 



i'" 1 A'-i 



(4.12) 



where jl\ = /ii(x, t). If jli were a function of time only, then we see that all the multipolar 
contributions in the right-hand side of ([4 .121) but the scalar / = one would be zero, because 
of the factor r'\ with I > 1 (this is clear with the old regularization, and easily checked to be 
true with the Lorentzian regularization as well). We defined Vdistr as being V but computed 
with the function of time [/ii]^ instead of the true /ii(x, t). This Vdistr is exactly the one 
which has been computed in the subsection p^V A| . [It can be checked that up to the 3PN 
order [p^i]-^ = (/ii)i-] Therefore, by the previous argument, Vdistr is produced entirely by the 
scalar part / = in the latter multipolar expansion, so that its complementary to the true 
potential V reads as 



V-Vr 



distr 



G 



+ 00 



\n Qn 



hoc 



^ n!c" df^ I ^ /! 

n=0 L 1=1 



+ 1^2 



(4.13) 



^ 1 



where the sum over / starts with / 

iliL- 



1. Thus, the problem reduces to the computation 
of each regularization terms Obviously, at a given post-Newtonian order, these 

terms will all become zero for I larger than a certain value. We find that, up to the 3PN 
order, all the regularizations are zero starting at / = 3, namely [r'ln'^jl'^^ = 0{8) for any 
/ > 3, while the non-zero values for / = 1, 2 are given by 



r in i/ii 



r 271 //il 



3 — n 



12 + 0{8) , 
m\m2 



c"ri2 



?,G'^m\ 



1^1 



Replacing these results back into ( [4.13|) , and using the fact that 
intermediate form 



(4.14a) 
C(8) . (4.14b) 
= 0, leads to the 



V -Vc 



distr 



12 



G^m\ 
2c6 



v\v{di. 



2c6 W 



C(8) + 1^2. 



(4.15) 



As we see, the non-distributivity of the potential \^ is a 3PN effect. Expanding the time 
derivative in the last term we find that the dependence on the velocity v\ cancels out, 
which is normal because a velocity-dependent term would violate the Lorentz invariance, 
in contradiction with our use of the Lorentzian regularization [F]^. The final expression is 
simple: 



V-V, 



distr 



5 G^m\m2 



n\2di 



12 



0(8) 



2 . 



(4.16) 
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The contribution of the non-distributivity in the acceleration of 1 is given by the gradient 

clIj X cLS 



W\ - [9.Kiistr]i = 5 " g \ ^ ni, + 0(8) . (4.17) 

^ ^12 

C. Computation of quadratic potentials {dVdV) 

By definition, the quadratic potentials are those whose sources are made of products of 
two compact factors, like V, Vi, W-j , etc. (or their derivatives, in most of the time). A 
typical source term for them is of the type ^^dVdV" , hence their denomination; for instance 



But the quadratic source terms may also involve other quantities of the same structure, as 
it is the case for dtW-^^dkV appearing in the source of the potential Y-^^^^'^'^ [cf ( |3.27b| )]. 



1. Matching to the external field 



The retardation of the compact potentials defining the metric of an isolated fluid can be 
expanded in powers of 1/c only in the "near zone" -Dnear of the source, at a distance much 
smaller than the typical wave length of the emitted radiation. The question then is how to 
incorporate in the post-Newtonian metric the no-incoming radiation conditions at past null 
infinity. We achieve this by performing a matching between the post-Newtonian expansion 
of the metric, adequate in the near zone, and its multipole expansion, valid in the region 
Dcxt exterior to the compact support of the source. Recall that for slowly moving sources, 
one can always choose -Dnear and Dcxt in such a way that their intersection is not empty: 
-Dncar H Z^cxt 7^ (scc e.g. fS^). The field h^'^ admits a multipole-type expansion Jidi^h^"), in 



the sense of [Q, at every spatial point x G -Dext- As a matter of fact, it is shown in that 
the multipole expansion of the exterior field (a vacuum solution of the field equations) that 
matches, according to the theory of matched asymptotic expansions, to the post-Newtonian 
expansion in the interior of the source, is given by 



-1 



B 



^4l_.'-^d,r-n^at-r/c)\ (4.19) 

1=0 



(with L a multi-index of order / G N). The multipole moments Ti^ entering the right-hand- 
side read as 



FPb^o jd'y (^M^ i/Lr^^(y,M) 
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where t^" represents the (formal) post-Newtonian expansion of the complete source term 



'T^'^ = Idl'^^'^ ~^ li^^^'^ of fi^l^i equations ( |3.3| ). These expressions are defined by analytic 
continuation in B, and the symbol FP^^o denotes the finite part when B goes to zero of the 
Laurent expansion of the analytic continuation (we refer to p^ , |5B| for more details about 
this finite part). 

Let us show how we find the "matched" solution of the equation OP = S at the rela- 
tive IPN order (this is all we shall need in this paper). We neglect all higher-order post- 
Newtonian contributions in the source term S, and look for the solution of 



□P = + C(3) . 



(4.20) 



Since the formula ( [4.19| ) results from the properties of the d'Alembertian operator (and is 
not specific to the field variable h'^'^), we can use it with the replacements of Ai{A) by 
M{S^™) and of r by j^a^^™- Thus, the multipole expansion of the solution must satisfy 



B 



MiS 



1PN\ 



1=0 



with 



;W(P) = FPB^o□7^' 
+ 0{3) , 



VL{t-r/c) 



(4.21) 



The partie-finie retarded integral of the multipole source A^(S'^^^) has to be handled with 
care. It is not licit to develop when c — > +oo the integrand under the integration symbol 
because the source is not compact supported. The correct formula was shown in |E6i to be: 



FP 



To J 



FP 



(2fc)! \cdt 



k=0 



2k 



dh 



-An 



X — X 



/|2fc-l 



B 



To' 



1>Q ^ 



nL{t-r/c)-TZL{t + r/c) 
2r 



(4.22) 



The hat on the partial derivatives Ol indicates that the trace has been removed, i.e. Ol = 
STF(9l). The IZl functions parametrize the general solution of d'Alembert equations that 
are smooth near the origin: "antisymmetric" solution as given by the last term in ( [4.221 ). 
We have, more precisely. 



with: T/(y,u 



TIM = FPs^o j d^'y (^^^ r;(y, u 
, ,,+i(2/ + l)!' 



(4.23) 



2'/! 



dz {z^-\)'M (5^™) (y,n-z|y|/c) 
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Here, denotes the symmetric trace- free tensor associated with y*^ . . . y*' , for / G N. With 
( [4.21| ) and ( ^.22| ) we can write the multipole expansion M.{P) at the IPN order as 



M{P)=YVb- 



dh 



-An X — x' 



B 



-] Mis 



IPN 



)--y 



1=0 



-9l (I ) V,{t) 



+ 



+ 



d 



2c2 V dt 
1 



FP 



B- 



Atic 



n{t)+v{t) +o(3) 



(4.24) 



where the last term, of order 1/c, is a simple function of time made of the functions 7^(t) and 
V{t) defined as being 7^L(t) and VLit) with / = (the dot indicates the time-derivative). 
Now, it can be shown that the latter multipole expansion can be re-written under the new 
form 



A^(P) = A^(p(i))+ n{t)+V{t) 



d_ 

2c2 



[^(P("))] +0(3) 



or, equivalently (indeed the second term is a mere function of time, and the multipole 
expansion obviously commutes with the time derivative). 



7W(P) = 7W PW + 



Atic 



In these equations, P^^-* and P^^^^ denote the matched solutions of the following Poisson 
equations 



^p(ii) = 2P(i) . 



(4.26a) 
(4.26b) 



Therefore, we have reduced the problem of finding the matched solution of the d'Alembertian 
equation ( [4.20| ) to that of solving and matching the two successive Poisson equations (f4.26aj) 
and ( ^4.26b| ). Now, from the equation ( [4.25| ), it is evident that the correct matched solution 
of (|4.20|) reads in terms of the matched solutions of (|4.26|) as 



p(i) + 



Atic - 



7^(^) + v{t) 



+ 



2c2 



[p{n)] + (9(3) . 



(4.27) 



To recall the meaning of this solution we shall often denote it as P 



for P = PjlJitch (fo^ instance fifmatch computed below) and P = Pj 
we shall find that the function TZ{t) appearing in the 1/c term of our solution ( [4. 271) is in fact 
always either zero or cancelled out by a spatial gradient in the case of the applications made 



(11) 

match 



-Pmatch below; similarly 
(e.g. /match)- Actually, 
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in the present paper. Thus, it will not be considered in this paper, whereas the function 
V{t) plays a role and is given by 



B 

r ' 



V{t) = Pf { FFb-,0 Id'^ {-] ^'™(^. t) } (4-28) 

(of course 5"^^^ there could be replaced with this approximation by 5"°-^^^) . See ( [4.41j ) below 
for an example of computation of this function. 

In practice, in order to find the matched solution of a Poisson equation, P^^^ for instance, 
we proceed as follows. Suppose that we know a particular solution of the equation, say Ppart- 
Then the correct solution is necessarily of the type P*^^^ = Ppart + ^hom' where h^l^ denotes 
an homogeneous solution of the Laplace equation (harmonic function): Ah^^^ = 0, which 
is moreover regular at the location of the source points. Note that its multipole expansion 
coincides with itself, A^(/ihom) ~ ^iom- Now, the latter homogeneous solution is determined 
by the matching equation as 

-M{pZ)- (4.29) 

It is not evident on that expression that the right-hand side is an harmonic function; but 
it really is, as can be verified explicitly in practice. We compute the multipole expansion 
of the source term S-^™ as well as of our particular solution Ppart- our case this means 
computing the formal expansions of S-^™ and Ppart when r tends to infinity or equivalently 
when the two source points yi^2 tend to zero. The computation is greatly simplified if one 
considers the dimensionality of the source. Suppose for instance that [5*^^^] = (length)*^ 
which means [h^J = (length) ^^+2. Then, using the fact that this function is harmonic, its 
structure is necessarily of the type h^^m ^ S ^ with I + li + I2 = d + 2. This shows 

that in order to obtain h^l^ completely it is sufficient to develop the right side of (|4.29|) 
when yi^2 — > up to the order d + 2 included (i.e. to control all the terms ^ in the 
expansions which have h + h < d + 2). All the higher-order terms, having li + I2 > d + 3 in 
the right side of ( [4.29|) must manage to give zero. The same method is used to compute the 



homogeneous solution h^^a contained in P^^^\ We shall implement this method in practice 



horn 

below. 



2. Structure of the quadratic sources 

In the context of the present paper, we will not need to compute the quadratic sources 
beyond the IPN order. As a consequence, we will deal with only a few kinds of elementary 
sources. By equation ( [4.7|) , we already know the structure of V at the IPN approximation. 
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The other compact retarded potentials have a very similar form. After expansion of the 
retardation of any of them, say P^'~'\ we get: 

Foi 



J — 47r|x — x'l c J —Atx 



+ -^dlVi [ ^|x - x'l F6i + 1^2 + C(3) , (4.30) 



2c2 * / -47r 



with F G J^. The first contribution has been calculated in [see equation (6.18) there]. 
What is interesting for us is that the result writes as a sum of space derivative of 1/ri 
(or I/T2), i.e. 9i(l/ri) (with the convention that L designates a multi-index of length /). 
Similarly, it is easy to convince oneself that the third contribution is composed of terms diri 
(or dLT2)- Moreover, the action of time derivatives in front of the integral leaves the latter 
structure unchanged, in accordance with formulas such as dtd^ri = —v\didLri. The second 
contribution in ( [4.30|) is a mere constant with respect to x. In fact, as they appear in the 



quadratic sources, the compact potentials are preceded by some space or time derivatives. 
Now, these derivations have to be performed in the sense of pseudo-functions . From these 



considerations, we are now in a position to tell what is the precise structure of the sources of 
quadratic potentials. They read as a sum of what we shall call elementary terms. As we are 
interested here in their spatial behaviour only, we shall omit purely time dependent factors, 
though they are normally included. Newtonian elementary terms are themselves products 
of two pseudo-function derivatives of contributions coming from the first integral in the 
generic expression ( [4.30| ): 9jPf9it(l/ri) x diPidMi^/Ti), or djVidK{'^/ri) x dLPidM{^/'r2) 
(and similarly with 1 <-> 2), where J, L, M are multi-indices of respective length j, 
k, I, m. In the same manner, the IPN terms result from products of pseudo-function 
derivatives of Newtonian and post-Newtonian integrals as the first and the third ones in 
( pq ): djVidK{l/ri) X dLPiduri, or 5jPf5i^(l/ri) x ^LPf^A/ra, and 1^2. As for the 
0.5PN terms, they are simply the pseudo-function derivatives of ^^(l/ri) or i9l(1/''"2) (times 
a mere function of time). It is also natural to distinguish between the "self" elementary 
terms on one side, which depend on one body only, e.g. di'P{{l/ri) x 9jPf(l/ri), and always 
admit pre-factors /i^, fiifii, and the "interaction" terms on the other side, involving both 
objects, e.g. (9jPf(l/ri) x 9jPf(l/r2). The 0.5PN terms 5i(l/ri 2) are considered separately. 
To be more explicit, we shall provide as an example the IPN source of W'i/''^'') defined 



by (|18D: 



^^(.vav) ^ _ G^~^l8,l>fl a,Pf 1 + ^ (af5,Pf 1 d,,,Fin - .^l^Pf ^ 5,)..Pfr,) 

- G^flifl2diPf— a-Pf — (4.31) 

ri r2 



Here, we have used the fact that PfdiVi = QjPfri and PfdijVi = dijPiri. 
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The sum of the retarded integral of the elementary terms then gives us the complete 
quadratic potentials after expansion in 1/c and matching. Therefore, these potentials are 
generated by the sources through some partie-finie integrals, which can be regarded as the 
result of the action of the elementary terms, considered as pseudo-functions, on smooth 
quantities in the field point. By inspection, it can be shown that the distributional part of 
the self terms never contributes to the previous integrals, whereas the partie-finie derivatives 
applied to the interaction terms coincide with those of the Schwartz distribution theory. 



3. Integration of the elementary sources 



We now come to the solving of DP = at the IPN order for each of the elementary 
terms composing the quadratic sources. We proceed following the method we exposed at 
the end of paragraph [IV C 1| . For this purpose, we first need to find a particular solution of 
the following Poisson equations: 



(I) 



and 



part 



ri r2 



AR 



(I) 



part 



1 

n 
1 

r2 



(4.32a) 
(4.32b) 



with L = ii . . .ii and K 



Ji ■ ■ -Jk- 



From Ppart) we deduce the matched value P^^^ by 



computing h^l,^ according to the relation ( [4.29| ) adapted for each elementary terms. 

The equation ( [4.32a| ) involves only the vector variable ri, so that it is simple enough 
to be integrable in a systematic way. To put the sources into a more suitable form, we 
start by applying the derivative operator that enters the self terms in the sense of functions, 
since the purely distributional part of the derivative does not contribute. The result is an 



adequate power of Vi times a finite sum of partial terms 5*^*^ . . . 5*^*:- 



1*2*: 



n 



*2fe+l 
1 



n 



1 5 



we shall 



denote more compactly as 6 



2K„L-2K 



n 



The solving of the Poisson equations rests then on 



the well-known identities (easily checked by direct calculation): 
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In ( ^ 1 + 



2/ + 1 
1 



for a G C \ {/ - 2, -/ - 3} 
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(4.33) 



where n\ is the trace-free part of raf, and rio a strictly positive constant. The quantities 
between braces are particular solutions -Ppai-t of AP*^^^ = r^nf, and we must in general add to 
them some harmonic functions to be evaluated by matching to the external field. Equations 
( [4.32b| ) are a priori the most difficult ones, because of the mixing of the sources 1 and 2. 
As a matter of fact, determining Pp^rt amounts to solving: 
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Ag = A/12 ^ n ^ ^4 34^ 

rir2 r2 



in the sense of distributions, on account of the fact that, for instance. 



ri r2 ' rir2 



where dn and d2K denote the partial derivatives with respect to yi and y2 [the same 
transformation apphes to diPiri x (9xPf(l/r2)]. As a consequence, lQk = <9iL<92xfl' and 
lIk = 9iLd2Kp'^ clearly verify: 

A l9k] = a^Pf-^^Pf- and A Lfk'] = ^xPfri^^Pf- • (4.35) 

'- -' ri r2 r2 

Note that the derivatives above should be understood as mere (Schwartz) distributional 
derivatives. Luckily, particular solutions of equations ( |4.34| ) in the whole space can be 
exhibited Ogg. We may take: 



g = \nS , with S = ri + r2 + , (4.36a) 

= -^riri2{nini2) (^g - ^ + ^('"2'"i2 + rir2 - riVu) , (4.36b) 

where (ninu) denotes the scalar product of Euclidean vectors. The function g is symmetric 
in its three variables x, yi and y2, so that 

1 1 

Aig = and A25' 



^^1^12 r-2ri2 



in the sense of distributions. For a more complete list of useful formulas, see |32|. We have 
also the two identities: 

AJ'' = 2g and A2/12 = ^ . 

r2 

Their proof is straightforward and call for some simple relations permitting to express the 
scalar products (nin2), (nini2) and (n2ni2) by means of some fractions involving ri, r2 and 
ri2. These relations, given by (5.14) in ||HT|, are very convenient in most of our computations. 
The function /^^ is obtained by exchanging x and yi in the function / which was introduced 



in the appendix B of ||57 



/ = ^rir2{nin2) (^g - ^ + ^('"I'^ia + ^2^12 - ^^^2) , (4.37) 
and which satisfies, in the sense of distributions, the equation 

A/ = 2^7 . 
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Once Poisson equations with IPN source are integrated, it remains to find the homo- 
geneous solutions to be added to get the full matched solution. Most of the self terms are 
already correct, namely those that go to zero when r — > +00. The other ones are determined 
from the interaction terms by taking the limit y2 — ^ yi, which happens to be always possi- 
ble. The matching formula ( |4.29| ) provides the function h^l^ associated to Ppart = 9- The 
computation is very easy because the dimension of the source is [l/(rir2)] = (length)"^ (i.e. 
d = —2 in the notation of the end of paragraph [IV C 1|) , therefore one needs to control only 



the constant term in the right side of ([4.29|) when yi_2 0. We arrive at h^l^ = — ln(2ro) — 1 
hence the correct fi'match solution of A^match = l/rirg: 

f S \ 

fi'match — In I 7; ~ 1 5 



2r, 



where Tq is the positive constant occurring in ( |4.19|) . Similarly, but with a little more work 



because the dimensionality of the source is now (i = so we must expand to second order 
in yi,2, we obtain the matched value corresponding to /^^ as 



/match = -^riri2{nini2) (^fimatch " g J + g (^^2^12 + ^irs - riris) - -r - - (yi^/a) + 2 ^(^2/2) 

(where (nyi) for instance denotes the scalar product of n = x/r with yi). As a consequence, 
the potentials P^^^ satisfying ( |4.32b|) are respectively (-)'+^Lfi'match k and (-)'"^^L/match k- 
With this result in hands, we are able to deduce very simply all the self terms that do 
not match properly yet. We shall content ourselves with examining how this works on an 
example. Let us suppose we want to solve: 



AP(^) = ndiiPf— . (4.38) 

ri 

We make the correspondence with the equation Af^^^^^ ^- = ri^yPf (l/r2), whose source 
coincides with ri^^Pf (1/ri) for yi = y2 (we recall our notation /^itch ij = ^2ii/mitch)- The 
distributional part of the derivative yields a compact supported contribution to /mLch ij 
given by 

< Pf n ( ^) . >= . 

\ 3 / — 47r|x — x'l 3r2 

which is zero in the limit y2 ^ yi, while the ordinary part yields a Poisson integral which is 
well-defined and easily evaluated for yi = y2. We conclude that the value of /^itch ij when 

y2 yi, 

P'^ — (/match tj)y2^yi ~ g*^ ^ ~ 2^l''^i ' 



is precisely the matched solution of the Poisson equation ( [4.38|) . 
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Let us complete now the program presented at the end of the paragraph [IV C 1| . Because 
of the presence of time derivatives at the IPN order, we restore in the elementary terms all 
the coefficients depending only on time, either through the trajectories yi^2 or the velocities 
vi 2 [we shall generically call a(t) this time-dependent coefficient; for instance a = (fif2)]. 
It is worth noting that the potentials assimilated to P^^^^ in ( [4.26b| ) are needed only at the 
Newtonian order, as they come with a factor. Consequently, the sources of the P'^^^^'s 
are simply (two times) the matched solutions of the Poisson equations AP^^^ = + 0{1), 
where are the Newtonian-type sources Q;(t) SiPf (1/ri 2) x S^Pf (1/ti,2); so, all we have 
to solve is: 

AP(") = 2a{t) rf6^'^n^-^^ , (4.39a) 
or AP^''^ = 2a{t) L9ra.tch K , (4.39b) 



with p G Z, and L, K some multi-indices. The elementary self potentials obeying ( |4.39a| ) are 



evaluated by application of the identities ( [4.33| ), before matching the full term (i.e. included 
the a coefficients) to the external field. The latter stage will be dropped here, because, 
on one hand, the general procedure has been explained before, and, on the other hand, 
powerful methods permitting to deal with the trickiest integrals one could encounter here 
will be expounded in Section 0. Let us look next at the second equation ( [4.39b|) . To get 
p(n) from some particular P^Ji one must perform the complete matching including all the 
time-dependent factors a{t). Here, for simplicity's sake, we give the result in the case where 
the (Newtonian) source is l/(rir2) hence P*^^^ = ^fmatch as we have seen before. Then, we 
have to find the matched solution of 

A particular solution Pp^rl of this equation is easily obtained with the help of the function / 
defined by ( |4.371 ) [indeed, the Laplacian of Ppart — / is a mere constant] . The corresponding 
homogeneous solution h^2n computed using the same equation as ( [4.29p but using the 
source 2P'^^) = Sgfmatch- The result, that we naturally call /match = P^^^\ reads as 

/match = ^rir2 (721^2) ^fi'match ~ ~ ^('^I'^is + r2ri2 - rir2) 

- \r{nyi) - \r{ny2) + ^{ym) ■ (4.40) 
o o 2 

Notice that in the case of the source l/(rir2) the only "odd" contribution 1/c in the formula 
(|]23) is that given by the function V{t) defined by ( [4.28|) ; the contribution due to TZif) is 
of higher order in this case. We readily find 



B 

;3, ' ^ 



V{t) = FPb^o d'^ [-] = -27rri2 (4.41) 

Jo) nrs 
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(no need of the symbol Pf). This calculation is also done in equation (5.8) of In Section 
|V|, we shall see more generally how such integrals can be obtained. Thus, our definitions of 
5'match and /match are such that 



( D^l^ ) = ^match - T^ri2 + ^T^^^ /match + C(3) . (4.42) 

\ ^^1^2/ match 

Finally, we have all the material to integrate the individual post-Newtonian terms in 
such a way that the inner metric matches to the external field at the 3PN order. Let us 
remark however that, in fact, the work we have done on the matching is, as seen a posteriori^ 
unnecessary. Indeed, summing up all the contributions in the potentials, we find that, had 
we made use of some "un- matched" elementary functions, e.g. g and /^^ defined by ([4.36|) , to 
compute the interaction terms instead of the corresponding matched quantities, and had we 
deduced jointly the corresponding self terms from the limit yi — > yi, we would have arrived 
at the same potentials up to the 3PN order. This means that the new contributions brought 
about by the matching to the external field actually cancel out in the final 3PN equations 
of motion. In particular, the constant tq which enters into the matched quantities (^match, 
/match /match disappears from the final result. Though we have verified this, we stick to 
our presentation and use systematically all the matched functions determined previously. 

To end this section, we shall achieve the example of the potential W^^^^"^ defined by 
( [4.18| ). We indeed already know its source from ( [4.31|) . We split the potential itself into: 



= -f^. - + - ■ (4.43) 



The first two contributions are respectively the matched solutions of the Poisson equations 

MJij = diUdjU and AKij = d(iUdj)d'fX , 

which come from diUdjU + ^d(^iUdj)d^X + 0{3) = diVdjV. Recall that the potentials U, X 
and so on have to be viewed as pseudo-functions (for instance U = PfG(/ii)j^/ri -|- 1 2), so 
the derivatives entering the source terms are distributional derivatives. The self terms can be 
determined with the help of the relations ( [4.33| ) and matching. To get the interaction part, 
we change the spatial derivatives to partial derivatives with respect to the source points yi^2, 
and next, we make the replacement l/(rir2) 5'match, '"i/'"2 ~^ /match' and rs/ri /match- 
The "odd" term Lij is a pure function of time given by 

= I —r- diUdjU , 





47r 


in 


il 



of Section |V|. The contribution Xij is the matched solution of the double-Poisson equation 



A'^Xij = 2diUdjU . 
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whose source is to be considered at the Newtonian order only. The iterative apphcation of 
( [4.33| ) plus matching yields the self terms; for interaction terms, we replace dL{l / rijdxi^ / ^2) 
by /match ii-- The results are 



' Sjlnri + — 



G'^jlljl2 iS'match 



] 1 



(4.44a) 



-^9,) Inr, + ^5^^5.1nri - ^9.,, (r^ln^) + ^ + §4 
A ■" 8 48 ^ ^ 16rf 8rf 



16 



H — -dii n ri + ^-^ 
16 ^ 96 



+ G'^^i/i2 [at fc(i/match j) + "^l"^! M(i/match j)] +1^2 



^?pf 



■Oi—di 



G^jlijl2dtdud2~ + 1^2 + 0(2) , 

^2,-2 



— 47r ri ri 



fiifi2diid2jPf 



dh 



1 



—Att rir2 



2 + 0(2) 



-dij (r^lnri) + d'^Hnri 



+ j/match J + 1 ^ 2 + const . 



(4.44b) 

(4.44c) 
(4.44d) 



In the last equation we do not write for simplicity a constant (associated with a function 
of type TV) which is cancelled out by the time derivative in front of that term. The self 
terms have been written in the form of some (ordinary) space derivatives in order to prepare 
the computation of the cubic sources. 



V. CUBIC POTENTIALS 

A. methodological scheme 

For methodological reasons, it is convenient to express all the cubic sources in a simi- 
lar way, with the help of the same set of elementary integrals. The so-called "cubic-non- 
compact" term 

^(CNC) ^ ^^^dvdv)Q^^y^ ^ ^5^^^ 

which is part of the X-potential [see ( |3.26a| )], is a good example to understand the successive 
transformation operations we perform in practice. Furthermore, this cubic-non-compact 
term is the only one we need to compute at the relative IPN order; all the other ones, which 
enter into T and 1^, are merely Newtonian. So the practical computation of ( ^.1|) is the 



most difficult one we face at the 3PN approximation. Recall that wj^^^^^^ was defined by 
( [4.18| ). We start from the expression of the source of X^'^^'^^ obtained by insertion of 
and d^) into (|3). We get: 



□^(CNC) _ _ ^ hi.dijU - ^KijdijU 
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^ [dlX,^d,,U + U,jd,jdlX] + 0(3) , (5.2) 



using the notation introduced in ( [4.43| ). In the right side, the potentials are seen as pseudo- 



functions (involving a Pf) and the derivatives are distributional. After carrying on the 
expansion of retardations up to the IPN approximation, we find: 



c^V 1 _ ^ _ I f 



y 47r |x — x'l c y 47r |x — x' 

^ ^'l^^l^ + - / ^ I n K 

c An / 47r X — x' 



+ ^a? y ^|x - x'l C7j,aj,[; + ia, y ^ Li,a.,[; + o(3) . (5.3) 

We have checked explicitly that the sum of the integrals occuring in this formula yields 
an integral convergent at infinity when considering the regularized value of the gradient 
{diX^'"^'"'') ^ which is the only thing required; thus we do not need to introduce a finite 
part at infinity (but of course the regularization Pf is needed to cure the point-particle 
singularities). The next step consists of replacing the potentials U, X, Lj^, Kij and Xjj 
given by (^4.44|) above by their values at the field point x 7^ yi and x 7^ y2. The spatial and 
time derivations appearing in each of the integrals of (|5.3|) are to be understood in the sense 
of pseudo- functions (see Section |I|). Consider, as an example, the term 



d^x 



' 1 



-47r|x-x'| 

Remind that Kij is given by (|]44^). Let us multiply (|4.44b|) by dijU = dijFf{Gfli/ri + 
Gjj'2/f2) and develop the product. The result is made of a sum of terms of the type (l/r^) x 
dijVi{l/ri), dk\nri x S^jPf (1/ri), dijki{rl\\ir2) x 9ijPf(l/ri), ikiif^^tcY,) j x 9ijPf(l/ri), etc. 
Some of them are functions of ri only; we call them "self" terms, whereas those depending 
on both ri and Y2 are called "interaction" terms. 



B. Self terms 

We agree on considering only the self terms (i) that are proportional to m\ rather than 
m\ (they are the same modulo the replacement 1 — 2), (ii) that do contribute to the IPN 
order at most. We leave aside the terms that are generated by Ljj, since their structure 
is especially simple and they are evaluated at the end of the section. By explicitly writing 
down all the sources, as done previously for KijdijU, we can draw the complete list of 
intervening terms. There are three types of terms: the VdVdV-tjpe concerns one kind 
of term only, i.e. 1/ri x c}jPf(l/ri) x djPf{l/ri); the so-called 3^-type refers to S^^Pfrf x 
dxPfrf terms, where p and q are positive or negative integers; the A/'-type terms come as 
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diPfdMiri In ri) x djPfdxrf (the terms y and Af are named after some integrals introduced 
below). There may exist contracted indices among the set of multi-indices. In particular, 
some terms involve a factor 5APf(l/ri) and are thus purely compact supported. In fact all 
the terms can be split into compact and non-compact parts. The latter part is an ordinary 
function that we are able to calculate explicitly. The former is determined from the results 
of Sections VI- VIII in and depends on the pseudo-function derivative we use. We shall 



refer to it as the self partie-finie-derivative contribution to the potentials. If we take the 
term dijkiirflnri) x dijPf{l/ri) for instance, it reads: 



1 1 

dijki{rl\nri) dijPi— = dijki{rllnri) dij— + dijki{rl\nri) Dy 



1 

In the case of the "particular" derivative defined by ( p.9|) , we have D?^'^*[l/ri] = 27rPf((5*-' — 
5n\n{)Si; so that 

d.Mrl Inri) d,,Pf- = 12 ^^^^, ^ + 87rPf ^^i^i/^^ (5.4) 

(like in the first term of the right side we sometimes do not write the Pf when there is no pos- 
sible confusion). In most of this section, we shall use the particular derivative D^'^'^*[F] given 



by ( p.9| ) instead of the more "correct" derivative Djj[F] defined by ( |2.1CI|) - (|2.11| ); in Section 



[VI| we shall discuss the effect on the final 3PN equations of motion of using the derivative 
Djj[F]. In order to obtain the self cubic potentials, all we have to do now is to apply the 
operator J |x — x'|^^ to the various sources we are focusing on, and J |x — x'| to 
the Newtonian source of X'^'-^'^^'-^). As a matter of fact, the resulting integrals can be viewed 
as partie finie pseudo-functions like (|5.4|) acting on l/(— 47r|x — x'|) or |x — x'|/(— 47r); both 



quantities are smooth at point 1, so the pseudo-functions associated with the non-compact 
part reduce to Schwartz distributions in that case (but, in order to construct the pseudo- 
functions themselves, we used the generalized distributions of 0). Each integral is indeed 
a sum of terms of the form r^5^^nf~^^, where p belongs to Z. It is convenient to write them 
as sums of pseudo-function (or, equivalently here, distributional) derivatives of quantities 
without indices ("scalars"), times some possible Kronecker symbols. We have for example: 

n{n{ _ 1 ^ 1 IS'^ 
rf 15 rf 5 rf 

= -d.jPi- + 7- + TF(8^>i - , (5.5) 

15 rf 5 rij" 15 r( 

and other similar formulas for n\/rf, n\n\n\/r\, etc. 

As an illustration of our handling of the sources, here are the effects of these trans- 
formations on dijki{r\\n.ri)dijPi{\/ri). Starting from equation ( |5.4| ) and using (p.5|), we 
find 



I ^:,M(r?lnr;) d'Pi\ 

J -An |x- x'l ^■'^'^ ^ r[ 
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where the first term is generated by the specific derivative ( |2.9| ). In this term the derivative 
can be changed to a partial derivative with respect to the point 1, and since we employ a 
pseudo-function derivative, we are allowed to permute integration and derivation symbols. 
This yields 



J -Att |x - x'I '^"'^ ^ r[ 



dh . 



5 J — 47r|x — x'I r'^ 5 J — 47r |x — x'| r^^ 



n-l'nl 36 5^' 



5r^^|x — x'I 5r^^|x — x 
The first two terms are left in this form for the time being. On the other hand the last term 



is computed following the procedure explained by the equations (6.17)-(6.18) in P? 



see 



also ( p. 16 ) below. By implementing the previous procedure for all the self terms entering 



KijdijU, we finally arrive at: 



Air |x 



30 ' J -47r|x-x'| 5 J -47r |x 

_ (n^ 7(W _ 1 ^ 2 + Oil) . (5.6) 

We follow the same way to treat the self parts of the other cubic potentials of interest here: 
^(CNC)^ y.(CNC)^ remaining of X(cnc). 

We find that there are definite contributions, coming at the 3PN order, due specifically 
to the pseudo-function derivative introduced in [^. Indeed, the distributional part of the 
derivative gives some well-defined non-zero contributions, while for instance the Schwartz 
derivative yields some terms which are ill-defined in this case. These contributions of the 
pseudo-function derivative actually take part in the values of Ti'^^'^^ and X^^^^^ only. De- 
noting them by ^seif^' and 5scif^ in the case of the particular derivative (|2.9| ) we find 

OseXiT = — — — (^1^12) + 1^2 (5.7a) 
iZ r^r^2 

, G^nv\ ( 17Gm2 , x 1 / x9 1 9\ -, „ x 

= (-^^'-■(".ni.) + 5;("i"0= - j + 1 " 2 (5.7b) 
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C. Interaction terms 



We consider exclusively the interaction terms (i) that are proportional to m\m2 rather 
than mim\, (ii) that contribute at relative IPN order, leaving aside those which are generated 
by Lij. Depending on whether they come from "simple" or "composite" cubic parts as shown 
respectively below, the elementary terms composing the sources read schematically 

9Pf[F(ri)] 9Pf[G'(ri)] dVi[H{v2)] (5.8a) 
and 9Pf [F(ri)] ^i^aPf [^(ri, ra)] = SaPf {9Pf [F(ri)] d{Pi[G{vi, v^)]} ; (5.8b) 

the functions G, H belong to JF, and it is also the case of F in general. However, there 
exist some composite terms for which F = Inri, but this is not a trouble since dF is still 
in JF. In the cases needed in this problem, G is always one of the four functions: fi'match, 
/match' /match' /match- Then, G is "regular enough" so that PfdG and dPfG coincide in 
any cases, and further simplifications of the sources do not seem to be possible at this level. 
All we need, thus, is to transform the simple cubic contributions ( |5.8a| ) similarly to the self 



terms (see paragraph |V lj| ). A typical example of elementary source we have to handle is 
l/(rir2) X dL^l/ri), where L represents at most two non-contracted indices. We can check 
that this term can be computed to Newtonian order only, hence it is given simply by a 
Poisson integral. In the language of pseudo-functions, this means that we have to evaluate: 

< -TT d'rPi^, —r^ r >=< 4 d'rPi^, , , rr > ■ 

r[r2 r[ — 47r|x — x'| r[ r[ — 47r|x — x'jrg 

The compact part of the dual bracket, which is associated with the distributional part of 
the derivative, when acting on l/(— 4'7r|x — x'|), i.e. 



<iD. 



1 



1 

> 



-47r|x — x'lrn 



{1=1 or 2), leads a priori to a non-zero result. Here, D^, denotes the distributional part 
of the multi-derivative, obtained in Section VIII of |47| and recalled by the equation (|2.12 ) 



above. However, the left side of the bracket, which is homogeneous to the (— / — 2)th power 
of a length, is necessarily of the type r|~'Pf5i, times some dimensionless angular function 
whose multipolarity differs from / by an even integer (because of the index structure of 
the operator D^^). Now, the previous compact part is equal to the angular integral of the 
r[~^ Taylor coefficient of l/(— 47r|x — x'|), times the angular dependence of Di[l/ri]. The 
integrand then appears as a sum of terms whose multipolarity differs from / + (/ — !) = 2/ — 1 
by an even integer and so, is always odd; thereby the angular integral gives zero. By similar 
arguments, we can prove that the other compact sources associated with the distributional 
derivatives will never contribute to the Poisson integrals constituting the potentials we are 
considering here. Actually, it is possible to put together all the various kinds of simple-type 
[see ( p.8a| )] cubic terms into a unique one, which is nf/rj"*"^ x 9Pf(l/r2). We express at last 



the first factor (if Z 7^ 0) as a sum of derivatives of "scalars" thanks to identities such as 
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( p.5|) . Since the pseudo- function derivatives will give here the same results as the Schwartz 
distributional ones, and by virtue of dLPi{l/rl)\v = Pf9L(l/rJ)|x) [where V is the set of 
smooth functions with compact support], the last transformation can be done in the sense 
of functions. Note however that the multiple derivatives of l/r2 are indeed distributional 
and play an important role in the sources. 

To sum up what precedes, all the interaction terms have the general structure: 



1 ^|x - xr d,LPf[d'jPfFir[) d,KPiG{r[, v',)] (5.9) 

(p = —1 or 1, F and G functions of JF). After commuting the integral and the derivative 
821, which is always licit for integrals converging at infinity (since 821 is followed by a Pf), 
the general cubic term ( p.9[ ) becomes: 

In the case where G = G{r2) (and k = 0), this rearranges as 

= {-ydij82LPi 1 ^|x - xr F{r[)G{r',) . (5.10) 
We shall end with implementing concretely our treatment of the source on two typical terms: 

/•dV^l 1 1 1 1;, 1 1 , 

/ —T] n — ^iPf- ^jPf- = o'9i*<92jPf / —-- , 5.11a 

J — 47r|x — x'l ■' r2 2, J — 47r|x — x'| r2 

/ -4vr |x-x'| ^^•^'"^t^h j ^ijPf^ = j |x-x'| ^'^^''^"^ ^'^^^V^ ' (5-llb) 



and by providing the complete interaction component corresponding to the Poisson integral 
of Kij8ijU, which completes the self part obtained previously in ( ^.6|) . The compact support 
terms have been explicitly determined, while the other ones are left un-evaluated for the 
moment: 

Wr) Pf / d^^' f21 n' 1 (^12«l) 

+a202ijfl j ^ fc /match '^ifc — 



^'■^ ' — 47r|x — x'l r^^r2 ' 16 J — 47r|x — x'| r'^ 
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96 



(i^x' 1 r'^ In r'^ 



dh 



—Att |x — x' 
1 



-An X — X' 



/I jklf] 



12 f,' 1 

match ^il / 



1 



—An |x — X' 
where D denotes the operator diid2i. 



7T '/i 



21 

match '-^il 



[ni2Viy 
8rf2r2 



1 ^2 + C(l) , (5.12) 



D. Elementary integrals 



1. Nomenclature 



The inspection of the formula ( ^.12|) for interaction terms issued from KijdijU suggests 
that we should re-express this potential, as well as all the other ones, by means of a restricted 
number of elementary integrals; basically one for each source type. Hence the proposal for 
a useful systematic nomenclature, which reflects their structure. We shall introduce the 
following notations (and ditto 1^2): 



y -- 

M 
1 



Pf 



Pf 



dh 



1 



—An |x — X' 



./I '1 '2 



^12 

1 (P,Q) 



dh 



1 



-Atx |x — x' 



S.pPf- 



7*^ ^' 1 'n 'Pl / ^-21 f 
■; T r-i mr, J- = 
-47r|x-x'| i(P,Q) J 



dh 



12 

match '^iP 



u = 

Ip 

K -- 
1 p 

Q 

1 (P,Q) 
1 (P,Q) 



47r |x 
rf^x' 1 



—An |x — X' 

rf^x' 1 



Tj- fcfl'match k 5pPf — 



—An |x 



7? fl'match C^pPf — 



s -- 



An |x - 
d^x' 



Pf 



— 47r |x 



> /match ^iP-Pf ^/ 



iQ/match 5jpPf — 



1 



—An |x — x'l 
rf^x^ 1 
—An |x — x' 
rf^x' 1 
— 47r |x — x' 



1 (">p) /" r/'^x' / 

.Q^7match ^ipPf ^ M=Vi j 3^|x-x'| ri^^Hn r'^ 



-An 
rf^x' 



X — X 



1 '2 



iQS'match 9pPf — 
'^l 

jQ^'match 



Q 

1 (P,Q) 



47r 
rf^x' 
-47r 



1 



X — X I jQ^fmatch f^jpPf " 



(5.13) 



The value of the previous integrals is not generally known at any space location, except in 
some special cases. The reason is that their sources involve three points, in addition to the 
integration variable x': the point x where the field is calculated, and the two source points yi, 
y2. The few ones that are computable for any x include notably some self integrals like y^'^'^^ 
and iS^'^'^-* , and the two integrals D and d2iQi (j,o) entering the interaction part of the 

XC-^^*-^) potential at the Newtonian order. Actually, there are no other cubic contributions up 



to the 2.5PN order, and that is why we were able in to get the complete expression of the 
metric in the near zone at this order. This property of the integrals D'^Mf' and d2iQi (i,o) 
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is linked to the specific form of the integrands, which are made of products of two second 
derivatives apphed on appropriate functions, such that the indices of the first derivative are 
contracted with those of the second one: Sj^Pf Inri x (5jjPf(l/r2) and {Qj x (9jjPf (1/ri). In 
both cases, particular solutions in the sense of distributions of the corresponding Poisson 
equations 



AKi = 2(9i,Pf Inri dijPf— , 



r2 



ri 



can be exhibited |]57| , |3T| . The solutions Ki and Hi of ( ^.14|) that go to zero as 



Ki 
Hi 



-A- Ai 



Inri 



r2 



+ 2 A, 



lnri2 
r2 



+ 



r2 



+ 



r\^r2 



(5.14) 
oo read: 
(5.15a) 



1 



-Ai 
2 

+ did2i 



q Inri 
— + 



n 



ri2 



ri + ri2 



-9 



lnri2 ^ Inri 



ri 



2r 



12 



d2i[idig)^] 

ri 



r2 



2r{r 



2^2 ' 



(5.15b) 



12 



with Ai = dm, A2 = d2ii] see also ( p. 321 ) for the expanded forms of these solutions. Thus, 
we have 

D^M'^'-'^ = ^Ki and ^2^gl ii,o) = . 



In the list ( p. 13 ) above, we note the appearance of iterated Poisson integrals such as 
A^(O'-i) and d2iQi (i,o) which come from the IPN contribution to the retardation expansion 
of xC^'"^^), and thus enter this potential through their second time derivative. What we shall 
have to compute for our purpose is their Hadamard regularized value. However, d2iQi (j,o) is 
not available on the whole space (for any x), so that we cannot deduce dfd2iQi (i,o) from it 
directly. We shall adopt then a different approach. In a first stage, we express the operator 
d"^ with the help of the partial derivatives du and d2i: if F(ri, r2) = d2iQi (4,0); then 



d^F = dt [v\duF] = v\v{dujF + ai^iiF + 1^2. 

Next, we commute the partial derivatives du with the integration sign, so that the deriva- 
tives act on the source of the Poisson integral. This operation is legitimate only if the 
derivatives diL when acting on the integrands are viewed as distributional (in the sense of 
Section IX in |^). The new integrands write then as a partie finie derivative of a prod- 
uct, e.g. diidgg d'ip^)] but remember that we are not a priori allowed to develop them 
according to the Leibniz rule in our formalism. In fact, these specific non Leibniz correc- 
tions happen to give zero contribution to the 3PN potentials. This can be seen by applying 
successively the formula (7.23) in [^] (which is indeed sufficient since the "test" function 
|x — x'l is smooth at the points 1 and 2) to all the sources we are dealing with. Therefore, 
we can employ the usual rule for derivatives of products to perform our final transformation. 
In summary, we will have for instance 
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/" (i^x' 1 f d^x.' 1 



(the index j means that j is excluded from the symmetrization operation). On this form 
we can apply the partie finie at 1 while staying in the same class of elementary integrals 
( ^■13| ). We conclude by going back to the example of the cubic term generated by KijdijU; 
we have, after appropriate reshaping, 

J -47r|x-x'| ^ ^ 1 6 30 ^ 5 



12rj 150rf A50rf 



I 4 1 48 1 1 (jj) ^1 



21 

(i,0) 



(-2,-1) .;2 (0,-1) (2,-1) 

Sriara 8 ^ 16 i 96 ^ i 

+^^^^^^^^^,^21 +i^2 + 0(l). 

i(i,jfc) i(i,o) 8riV2 J 



We refer to for the expressions of the other non-linear potentials expressed in this manner 
by means of the same nomenclature. The problem is now to evaluate all the elementary 
integrals from which the 3PN cubic potentials have been built. 



2. Parties finies of the elementary integrals 

As mentioned before, in most cases, we do not have at our disposal the explicit values 
of the elementary integrals in all space. This does not matter since all we need is their 
Hadamard partie finie at point 1 (or 2). Notice that the partial derivative with respect to 
y2 is the only one that commutes with the partie finie operation at 1; to be more explicit: 

[52iF(x,yi,y2)]i = 92i[F(x,yi,y2)]i , 
but [9HF(x,yi,y2)]i ^ 9h[F(x, yi, y2)]i (with F, G G ^). 

Thus, for each elementary integral, we shall determine first the partie finie of the quantity 
figuring under the derivation symbol 82 and, only then, apply the latter operator. On the 
contrary, we cannot bring the derivatives with respect to yi out of the partie finie at 1, so we 
are led to incorporate them to the sources, in the sense of pseudo-functions, by permutation 
with the integration sign. As a consequence, the integrals we are interested in are of the 
type: 

y 3^|x-xr 9i[PfF(ri)] 9i[PfG(ri,r2)] 
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= -—< 9i[PfF(ri)] ai[PfG(ri, r,)] , |x - x'^ > . 

They involve both a compact part (C) and a non-compact part (NC). The compact part is 
produced by the purely distributional contributions of derivatives in the integrands: 

< Di[F(ri)] PfaiG'(ri,r2) , jx-x'^ > < Di[G'(ri, rs)] PfdiF{ri) , Ix-x'^ > . 

47r 47r 

As mentioned at the beginning of the paragraph [V Q , the partie finie derivatives reduce 
here, in our case, to those of the Schwartz theory. This is obvious when G = G{r2), 
because then the source F{ri), regarded as a linear functional PfF(ri), acts on a function 
which is smooth in a neighbourhood of x = yi; in the other cases, the result follows from 
explicit calculations. Note that the number I of derivatives in front of G(ri,r2) = g'match, 
/match) /match) OT /match IS always Small cuough SO that 9iL_iG'(ri,r2) is bounded, hence 
(9iLPfG'(ri, = Pf(9i2,G'(ri, r2). Once we have in hands the compact part, it remains to 
obtain its regularized value at point 1. As a matter of fact, if the source is of the type PfF52, 
where F G JF, then 



< pf(f.,)(x') . ^ >= -i^^Ll! (.;^„-f) a, (1) (5,16) 

' ' l>0 ^ 

is smooth at point 1 and we need not call for the Hadamard regularization. Therefore, we 
are allowed to replace directly | x — x' | by in the left-hand-side of (|5.16| ) . When the source 
is of the generic type PfF5i, the same identity as ( |5.16| ) holds, but with 2 — 1. This shows 
that the integral < Pf(F5i) , |x — x'|^^ > is purely singular as ri — » 0, which means that 
it has no partie finie at x = yi: the 6i type sources do not contribute to the potentials 
computed at body 1. Summarizing, we have 

< Pf(F52) , -^-^ >] =< Pf(^'^2) , - > , (5.17a) 
|x — x'l ri 

< Pf(F5i) , = . (5.17b) 
1^ ^1 /I 



We refer the reader to equations (6.17)-(6.20) in |^ for more details. 

Let us now focus our attention on the non-compact parts of the elementary integrals, 
whose integrands are made, by definition, of ordinary functions. The problem is to get the 
Hadamard partie finie of the Poisson integral P oi F & T: 

PU) = Pf / FM , 



as well as the one of the iterated Poisson integral Q, given by: 

g(x') = Pf / ^|x-x'|F(x). 



45 



[Each of these functions depends also on the source points yi_2 and velocities vi^2-] Now, 
the partie finie prescription applies only to functions admitting power-like expansions near 
their singularities [see ( |2.1| )], whereas P or Q may contain logarithmic coefficients in their 
development: if we take for instance F = 1/rf, we shall have P = —[1 + ln{r[/ si)]/r[, where 
Si is the constant appearing in ( |2.3| ). Following we shall simply include the possible 
logarithms (i.e. Inr'^^) appearing in the zeroth power coefficient of the expansion of P or 
Q in the definition of the partie finie: see the equation (5.4) in [^. With this generalized 
notion of "partie finie", we can give a sense to (P)^ and {Q)^, as well as their gradients 



(diP)^ and (diQ)^. We make then the following statements (see Section V of [0] for proofs 
and discussion): 
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(5.18a) 
(5.18b) 
(5.18c) 
(5.18d) 



where the first terms in the right sides are made of some partie-finie integrals in the sense 
of the definition ( p.3| ). The "constant" r[ is the variable which tends toward zero when 
evaluating the partie finie. It is easy to show that the constant si cancels out between the 
two terms in each of the second members of the identities ( ^.18 ). Indeed, using the general 



dependence of the partie-finie integral on the constants Si, S2 as given by (4.20) in [^ 



we easily see that {P)i given by (|5.18a]) , for instance, depends on the constants r[ and S2 
through the formula 



(P)^ = -ln(^ (r?F),-ln +... (5.19) 



S2 I \ri 



2 



where the dots indicate the terms that are independent of any constants. As we see, the 
constant Si has been so to speak "replaced" by r[. This makes clearer why it is convenient to 
keep the In r[ in the definition of the Hadamard partie finie; if we had decided to exclude this 
logarithm from it, we would have found some bare lnri2 in the first term of ( 5.19|) instead 
of a nicer logarithm of a dimensionless quantity: ln(ri2/r'j^); but this is simply a matter of 
convenience, because we shall see that in fact the "constants" Inr'^ and lnr2 can be gauged 
away from the 3PN equations of motion. The same argument is valid for all cases in ( ^.18|) . 
As a consequence, the acceleration of the first body will depend only on two unspecified 
constants: Inr'^, and of course lns2 (and ditto for the acceleration of the second body). See 
Section |V11| for further discussion of these constants. 

The relations ( p. 181 ) answer the problem of evaluating the elementary integrals at the 
location of particle 1 without knowing their values at an arbitrary field point. The subsequent 
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task consists of calculating the partie-finie integrals in the right-hand-sides, which will turn 
out to be always possible. 



3. Integration methods 

The non compact parts of the regularized elementary integrals consist essentially of some 
integrals Pf J (i'^x F{ri, T2), with F & J-'. It is worth noting that the sources depend on ri, 
r2 exclusively, and not on the separate variables x, yi, and y2, because scalar products like 
(xyi), {xy2) or (^1^/2) occurring in f^^^^^,, /match' and /match are killed by the derivatives 
that precede them in the integrand. In this work, we make use of two different integration 
methods: (1) the angular method, in which we determine successively angular and radial 
integrals in spherical coordinates, (2) the analytic continuation method, based on the so- 
called Riesz formula. 

a. Angular method. Let F = F(ri,r2) be a function in the class JF. We assume for a 
moment that F is locally integrable at point 2; so we are allowed to compute the integral 
over the whole space, deprived from a small ball Bi{s) of center yi and radius s > 0. Thus, 
we start with the well-defined quantity /iR3\gj(5)'^^x F(ri,r2), but for convenience we write 
it in spherical coordinates (ri,^i,0i), such that the azimuthal angle 61 coincides with the 
separation angle between yi2 and ri, 

/ d^^ F(ri,r2) = / dn rl / dUi F(ri,r2) . 

Jm?\Bi{s) Jri J 

Actually the function F may be a tensor with many indices, but the only unit vectors in the 
problem are rii, 112, ni2, and only two of them, say rii and ni2, are independent, by virtue 
of the relation riiii + ri2ni2 = r2n2; and therefore, F can be expressed under the form of 
a finite sum of tensorial products of type nf'nfl (/, G N); moreover, each factor admits a 
symmetric trace-free (STF) decomposition on the basis raf = STFnf. Hence we have: 

F = J2J2^i''i2^m{n,r2) , lo, heN. 

1=0 k=0 

Here, the G'(ife)'s are some scalar functions of ri and r2. Now, r2 is related to ri and the 
scalar product (niriu) by r2 = \/rl + + 2riri2(nini2). So that the angular integral of F 



can be obtained by means of the "mean formula"; see, e.g., the formula (A26) in [Q. We 
get 



/ 'o feo pi / I \ 

F = 27r^^rif2<2 / dz G(ik) (n, y^r? + r?^ + 2riri22 j Pi{z) , (5.20) 

/=0 fc=0 "^"1 ^ ^ 

where Pi is the Legendre polynomial of order /. In the cases of interest here, G(ik) is always 
a sum of rational fractions with general structure r^rl/iri + r2 + ri2)", with p, q positive or 
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negative integers and n G N . It is easy to check that the result of the angular integration 
depends on the relative positions of ri and ri2. Therefore we must split the radial integral 
into two contributions according to the integration domains ri G ]0,ri2] or ri G ]ri2,+oo[. 
Typical terms coming from the angular integration are / (ri + ri2)", as well as some more 
complicated logarithmic terms such as r^ln(l + ri2/ri). Most of the corresponding radial 
integrals are obtained straightforwardly using some integrations by part, applying the partie 
finie at the bound ri = s [i.e., removing the poles with k > 1, and replacing Ins by 

In Si]. In the case of the latter logarithmic terms with p = —1, integrating by parts does not 
lead to anything, but the radial integrals can be found in standard mathematical textbooks: 
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(5.21) 



It can be shown that the integrals diverging at yi involve in general a logarithm ln(ri2/si) 
but never any vr^ terms. The procedure we have just described indeed permits calculating 
most of the elementary integrals. Consider for instance the integral 



1 (Jfi) 
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(5.22) 



(in which we replaced (^match by (7 = InS" since they merely differ by a constant). We are 



interested in the value of this integral at point 1, following the regularization. From (|5.17b|) 
we know that the distributional part of the derivative will not contribute. Then, using 
( ^.18a| ), we readily find 
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The non-compact integral has a sole divergence at point 1 , so that we can apply the previous 
method without any change, and get 



Pf 



c/^x 1 ^1 

i9 dj — 

—Air ri ri 



12"'12 



18ri2 12ri2 Gr^ 
As expected the constant Insi cancels out and we arrive at: 
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(5.24) 



A few elementary integrals diverge at the locations yi and y2 of both particles. In this 
occurrence, the integral Jj^a^^^^^^c/^xF has no meaning, and the previous integration process 
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is not adequate anymore. However, Proposition 2 in Section IV B of allows us to extend 
it to this case. We introduce the auxiliary source 



lb , 

2 



F2 = F-Y^ Ah 

6+3<0 

which is locally integrable near y2 but does not converge at infinity. As before, the angular 
integral of F2 around yi takes a different expression depending on whether ri < ri2 or 
fi > 1^12, so we must split the radial integration into the two domains ]0, ri2[ and ]ri2, +oo[. 
Then, with full generality, the partie-finie integral of the source F is given by 



rri2 r ^ r+00 

Ffs,,s2 I (f^F = Pf,, / dri rl / d^i F2 + dr^ rl / d^i 

J Jri2 



' 1 



+ 4.HF),^l„(?ii) 



(5.25) 



If this integral comes from a Poisson integral evaluated at 1, the constant si will be cancelled 
out as we have seen previously and replaced by r'^\ but there will remain in general a constant 
S2 coming from the singularity at the other point. With the angular method we were able 
to obtain all the elementary integrals (and their gradients) at point 1. See an appendix of 



54 for the complete list of those results. 

h. Analytical continuation method. The equivalence between the Hadamard partie finie 
prescription for integrals and the analytic continuation regularization has long been known 
(see e.g. [0), and we have recovered it in the three-dimensional case by the Theorem 2 of 



More precisely, for any F ^ T that behaves like o(l/r^) when r +00, the integral 
Jti'^x [txI S\)°'{r2l S2)^F of two complex variables a and /? admits an analytic continuation 
in the neighbourhood of a = /? = 0, and we have 
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(5.26) 



where FP means taking the finite part of the Laurent expansion of the (analytic contin- 
nation of the) integral when a — and /? — successively. This result is particularly useful 
in the case where F is of the type rfrg, with p and q relative integers, since the integral is 



directly computable thanks to the Riesz formula 51 : 



/■p / a+p+3 \ -p //3+2+3\ -p / a+/3+p+<?+3 \ 
axr^ r2 - vr ^ ^ ^ r^2 ■ \^-^') 

One may consult |^ for an example of practical computation. Most of the time, the 
structure of the sources of elementary integrals is more complicated than a simple r^r|, 
notably it involves many "free" tensorial indices which imply that generally the sources 
of elementary integrals, when fully developped, involve numerous inverse powers of S* = 
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Ti + r2 + ri2. However, by considering all the possible contractions of these free indices 
with vectors n\2 and Kronecker symbols (5*-', it happens that we reduce the computation to 
that of several scalar integrals we can obtain thanks to the Riesz formula ( p. 27] ) (i.e., when 
performing the contractions and after simplification of the result, we are always led to the 
simple structure r^'rf without 1/5 powers). It should be noted that the set of scalar functions 
that we compute contains the complete information about the complicated tensorial integral, 
i.e. it permits to reconstitute it exactly. 



Let us illustrate the method with the computation of the integral ( |5.23| ). It involves two 
free indices ij and is necessarily of the type 

/ ^ 7^^3d~ = <P{ru)n\2n{2 + ^{ru)S'^ , (5.28) 

where and ip are some unknown "scalars" depending on ri2. By contracting successively 
the integrand with n\2n{2 and 6"^^ , and simplifying, we find: 



gd — = — + — !^ - ^2 

12 12 « J Arirf2 4r^rf2 4rfr^2 



+ ^2 ^ ^ ^ 

4r^r^2 4r^ri2 4r^ri2 4r^ ' 
m ^ ^1 ^ ^2 



thus obtaining a sum of terms of the type r\rl which can all be integrated with the help of 
the formula ( p.27| ). This computation yields a system of equations for the scalars (p and ip: 
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By solving the previous system, and inserting the results into ( |5.28| ), we recover exactly 
the value given by (|5.2ij|). We have checked in this manner all the elementary integrals 



previously computed with the angular method. 



4- Finite part of integrals diverging at infinity 

For the moment, we have left aside the case where our integrals are defined by means 
of a finite part dealing with a divergence occuring at infinity. Such a study is needed to 
compute some 0.5PN integrals encountered in Section |IV|. From now on, we suppose that 
the source F admits an expansion when r — ^ +oo which is made of simple powers r"~^, with 
n < nniax (so the integral a priori diverges at infinity when nmax > 0), and we consider the 
quantity: 
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(5.29) 



We split the integral above into two integrals /jnt and /ext extending respectively over a 
domain Dint including the two local singularities yi,2, and the complementary domain Dext 
comprising the regions at infinity. Using the integration variable ri = x — yi, the external 
integral (on which the regularization Pf can be removed) reads as 

/e.t = FPs^o / rf'ri ("^±±1^^ F(ri + yi) . 

If we assume that the original integral can generate only simple poles ~ 1/B at infinity 
(which will be always the case here), we can replace it by 




+ ^^^-°fil "^'"'fe) ^^(^ + 2(niyi)^ + f|)i^ (5.30) 

Indeed, in the case of simple poles, the other terms, involving at least a factor B^, will always 
give zero. The second term is non-zero only if the corresponding integral does admit a pole, 
and can be calculated in a simple way by picking up the term of order 1/rf in the expansion 
of the integrand when r +oo. Notably, in the important case where the function F 
goes to zero like at infinity, its product with the log-term behaves at least like and 
therefore the second term in ( |5.30D gives no contribution. For instance, an integral divergent 
at infinity that we encounter in the problem is 



E. Lorentzian regularization of potentials 

All the potentials and their gradients (compact-support, quadratic and non-compact 
potentials) which have been computed in this section and the previous one were obtained at 
the point 1 using the standard Hadamard regularization {F)-^^. However, this regularization, 
being defined within the hypersurface t =const of the harmonic coordinates, must break 
at some point the Lorentz-invariance properties of the potentials. That is, if a potential 
defined for a smooth "fiuid" behaves in a certain way under a Lorentz transformation, we 
expect that its regularized value at the point 1 in the sense of (F)-^ will generically not 
behave in the same way. Nevertheless, the equations of motion in harmonic coordinates, as 
computed with the regularization {F)-^^, are known to be Lorentz invariant up to the 2.5PN 
order |3l|. Perhaps not surprisingly because of this fact, it turns out that the Lorentzian 
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regularization [F]-^ (defined in [Q) yields no difference with respect to the old regularization 
(F)^ for all the 3PN potentials but for one, namely the cubic-non-compact potential X*^^'^^^ 
defined by ( |5.1| ) and which had to be computed at the relative IPN order. (Evidently we 
have [F]-^ = (F)-^ for all the potentials which are to be computed with Newtonian accuracy.) 
Thus all the results obtained so far but the one for {diX^'^^'^'')^ [the relevant quantity for 
the equations of motion] are valid in the case of the new regularization. In this subsection 
we compute the remaining part [9,X(CNC)]^_(5^x(CNC))^. 

Since the regularization [F]-^ brings some new terms with respect to the old one [F)-^ 
starting at the relative IPN order, and since X^^'"^^^ is to be computed at the IPN order 
only, it is sufficient for this calculation to use the lowest-order, Newtonian value of X^'^^'^\ 
From the computation in ||57|j3T[1 we know the analytic closed-form expression 

of X(CNC) at 

Newtonian order for any field point (t,x). 



2 

171^1712 



1 



8r2rf2 16 



Oil) 



(5.31) 



where the functions Ki and Hi, which are solutions of certain Poisson equations, are given 
explicitly by ( ^.15 ). The complete developed forms of these functions are 
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We replace these expressions into ( ^.31 ), and we implement all the rules for the new regular- 
ization [F]^ defined in Section III of Equivalently, since the order of the computation 
is limited to IPN, we can use the closed form formula 



[F], - {F\ = - (n.v 



dtF+-v\d,F 



+ 0(4) 



(5.33) 



derived in Section IV of 1 48 1. As a result, we obtain, for the potential itself. 



G^m?m2 (4:3, ,o , , , , 43 n 1, , , , 
= ^24 |i^(^i2t;i)' - {ni2Vi){nuV2) - —vl + -{V1V2) \ + 0{4) . (5.34) 

In the case of the gradient needed for the equations of motion, we get 
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Y^{ni2Vi)vl + —i7ii2V2)v{ + — (ni2t;i)t;^ [> + 0(4) 



(5.35) 
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As we see, the new regularization brings some definite non-zero contributions at the IPN 
order in the case of this potential, which will constitute a crucial contribution to the 3PN 
equations of motion. The right-hand side of ( |5.35| ) is not invariant by itself under a Lorentz 
transformation — it cannot be — , but will ensure finally the Lorentz-invariance of the 3PN 
equations of motion. 



VI. LEIBNIZ TERMS AND NON-DISTRIBUTIVITY 

A. EfFect of a gauge transformation 

In this subsection we study the effect of a gauge transformation on the 3PN equations 
of motion as well as energy of the two particles. Let {x'^} denote the harmonic coordinate 
system and g^u{x) be the harmonic-coordinate metric, generated by the two particles, that 
we have iterated in previous sections up to the 3PN order. The metric depends on the 
position X of the field point, and on the coordinate time t = x^/c through the trajectories 
yi,2(^) and velocities vi^2(^) of the particles, i.e. 



g^^Ax) = 5'Mi/[x;yi(t),y2(t); vi(t),v2(t)] . (6.1) 

We know that the dependence of the metric over the velocities arises at the IPN order 
(see e.g. the equations (7.2) in [0), namely the order (9(4,3,4), where this notation is 
a shorthand for saying 0(4) = 0{l/c^) in goo, 0{3) in goi and C(4) in gij. Consider an 
infinitesimal coordinate transformation of the type 

= + ^^'{x) , (6.2a) 
e(a;)=e[x;yi(t),y2(t)], (6.2b) 

where, in order to simplify the presentation, we assume that the gauge vector depends on 
the positions yi,2 of the particles, but not on their velocities. Furthermore, we suppose that 
this gauge transformation is at the level of the 3PN order, which means that = 0{7) and 
= 0{6), or equivalently = (9(7,6). In addition, in a first stage, we suppose that the 
vector ^'^(x) is a smooth function of the coordinates even at the positions of the particles. 
The new metric in the new coordinate system {x'^} is 



= 9'^,A^'■,y'lit'),y2it')■,v[{t'),V2{t')] , (6.3) 

where the new trajectories and velocities y'^^ 2) ^12 parametrized by the new coordinate 
time t' = x'^/c. The coordinate change (6^), when applied at the location of each of the 
particles, yields the relations between the new and old trajectories, which, when retaining 
only the terms up to the order 0{6), read as 
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y'\{t')=y\{t)+e{yi) + 0{8), 
y'l(t')=ylit)+C{y2) + Oi8), 



(6.4a) 
(6.4b) 



where the ^^(2/1,2) 's denote the gauge vector at the position of the particles, for instance 
^^{yi) = ^'^[yi(^); 71(^)5 y2(^)]- The new metric ( |6.3| ), when expressed in terms of the old 
variables, follows from this as 



g'^ui^') = %u{x) + C{x)^^g,. + C{yi)Ag,u + C{y2)2d,g,. + 0(10, 9, 10) , (6.5) 

where we have used the fact that the dependence of the metric on the velocities starts at 
the IPN order, so the terms due to the modification of the velocities do not contribute to 
( |6.5|) . Since the 3PN metric depends on space x only through the two distances x — yi and 
X — y2, we have dig^j_y + idig^u + 2dig^u = 0, and so an equivalent form of ( |675D is 

gU^') = V(a^) + ICiyi) - C{x)]l^^g,u + [e{y2) - e{x)]2d^g,u + 0(10, 9, 10) . (6.6) 

The equation (|6.6| ), when combined with the law of transformation of tensors, i.e. in the 
present case 



g,Ax) = g'^^x') + d,^, + d,^, + 0(10, 9, 8) , (6.7) 

where = rj^^^", gives the metric variation or Lie derivative S^g^^, = g'^^{x) — g^iv^x) (where 
the same variable x is used for both the transformed and original metrics) as 



hg^^u = -d,^u - d,i, + [C{x) - C{yl)]l^^g,u + [C{x) - C{y2)]2dig^.. + 0(10, 9, 8) . (6.8) 

In fact, up to this order, only the 00 component of the metric includes a "non-linear" 
correction term; and, within that non-linear term, the metric can be approximated by its 
Newtonian part, so 

5^^oo = -290^0 + '^{\C{x) - C{yi)]AU + [C{x) - C{y2)]2^^U^ + 0(10) , (6.9a) 

kg^i = -do^r - 5.^0 + 0(9) , (6.9b) 
hg^j = - + ^(8) , (6-9c) 

where U = + is the Newtonian potential (with a small inconsistency of notation 



with respect to previous sections). Now, it is easy to check that, in the sense of distributions, 

^(iCix) - e{yi)]idiU + iCix) - e{y2)]2diU^ = -2diQd,,u - a^au . 
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Indeed, the delta-functions at the points 1 and 2, which come from the Laplacian of [/, are 
killed respectively by the factors — CiVi) and — C{y2), which vanish respectively 
at these two points, in front of them. So, we can write for 5^(?oo the simpler but equivalent 
expression 



%oo = -290^0 - ^A-'[2d,^^di,U + A^U] + 0(10) , (6.10) 



where denotes the usual Poisson integral. 

The latter result can be generalized to our framework of singular metrics by allowing 
the gauge vector to become singular at the positions of the particles (in the sense that 
^'^ G JF), provided that the integral appearing in ( p.lQ ) is treated as the Hadamard partie 



finie of a Poisson integral in the way which is investigated in Section V of Let us 

consider, for example, the 3PN gauge vector given by 



^/^ = — — + — 



(6.11) 



where ei and €2 denote two dimensionless constants or possibly functions of time t (and 
where m = mi + 1712). Note that with this choice of gauge vector the new coordinates 
satisfy the condition of harmonic coordinates outside the singularities (i.e. in the sense of 
functions) at the 3PN order: indeed D.^^ = (9(9,8). When inserting ( |6.11| ) into ( |6.10| ), we 
must be careful about evaluating the last term of ( |6.10| ) in the sense of distributions, taking 
into account the fact that A^j is distributional. For this term we obtain 



nj \r2 



where 7^ and 72 are the Newtonian accelerations of 1 and 2. Therefore, we find 



2G^m^ 



+ 2A-^ 



dijU 



+ 0(10) + 1^2. (6.12) 



In the case where ei and €2 are some pure constants (independent on time) we can somewhat 
simplify the latter expression by using the fact that the accelerations cancel out in the first 
term. In this case, we obtain the full metric transformations as 
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2G^m^ 
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ri 
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2G^m^ 



,7 '''Aj ( ^ ) + ^(9) + 1^2, 



c' 



,6 



0(8) + 1^2 



(6.13a) 

(6.13b) 
(6.13c) 
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By comparing this with the 3PN metric ( p.24| ) , we see that the gauge transformation induces 
the following changes in the 3PN potentials T, Yi and Zij: 



5,f = - 






16 




G^m^ 
8 




G^m^ 

Q 


5^ Zij - 



-1 



+ 1 



1^2. 



(6.14a) 

(6.14b) 

(6.14c) 



The computation of the non-linearity term in ( |6.14aD is straightforward, and we get 



dijU 



Grrii 
'2^ 



Gm 



2 ij9ij 



where g = ln(ri + r2 + ri2) is a kernel satisfying Ag = [see also ( [4.36| )], and where we 



denote ijgij = diijd2ijg = D^g. At last, we insert the latter changes of the 3PN potentials 
into the (regularized) equations of motion p.32| ) with ( p.35| ) (there is no need to include 
a correction due to the non-distributivity), and obtain the corresponding change in the 
acceleration of the particle 1 as 



2G^m^ 



12 



617712 - e2mi n 



''12 



+ 



G^m^e2 



12 



- 15(ni2fi2)^<2 + 3f i2<2 + 6(ni2t;i2)t^ 



12 



(6.15) 



In the case where ei and €2 depend on time, there are some extra contributions proportional 
to 62 and 62. A good check of ( |6.15| ) is the fact that to the change in the acceleration (|6.15|) 
always corresponds a change in the associated energy; that is, the gauge transformation 



does not modify the existence of a conserved energy (see Section [VI 1| for the computation 
of the 3PN energy). Namely, we find that the combination mi5^a\v\ + m2S^a\vl is a total 
time-derivative, and from this we obtain the gauge transformation of the energy as 



G^ 



£2- 



c^r 



12 



Gm2 
ri2 



3(ni2fi)(ni2fi2) + (fifi2) 



+ 1^2 



(6.16) 



B. Leibniz contributions 



An important ingredient of the present computation is the novel distributional derivative 
associated with the Hadamard regular izat ion which has been introduced in (see also 
Section This derivative permits us to derive in a systematic and consistent way all 
the integrals encountered in the problem; however it represents merely a mathematical tool, 
which is maybe not connected to any relevant Physics. Therefore, it is important to know 
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exactly the role played by this derivative in the 3PN equations of motion, with respect to 
say the Schwartz distributional derivative We know that our distributional derivative 
affects the computation of two types of terms: (i) the "self" terms entering a priori in the 
non-linear potentials X, T and Yi and which are ill-defined in the case of the Schwartz 
derivative (see Section |V]), (ii) the "Leibniz" terms which account for the violation of the 



Leibniz rule during the 3PN iteration of the metric as discussed in Section In the 
present subsection, we compute the Leibniz terms and combine the result with the one of 
Section concerning the self terms. The conclusion is that the terms coming from the 
use of the distributional derivative are necessary for keeping track of the Lorentz invariance 
of the equations of motion, and that no other Physics is involved with them in the present 



formalism (see also Section [Vll| ). We do the computation for both the "particular" derivative 



defined by ( |2.9| ) and the more correct one given by ( p.lO| )-( pTTl| ). 

The Leibniz terms discussed in Section |ITT| consist of those contributions of the type 
( 3.23| ) and alike which arise in the process of simplification of the 3PN metric h^" by means 



of the Leibniz rule. These terms depend only on the distributional part of the derivative, 
Df'*[F] or Di[F]. The formulas giving the complete Leibniz terms in h^^ , not being very 
attractive, are relegated to Appendix 0. When reducing explicitly these formulas we find 
that all the terms take the same simple structure and, not surprisingly, arise only at the 3PN 
order. As already announced in (p.27|) , the Leibniz terms imply a priori a net contribution 
to the 3PN potentials T, Yi and Zij. Actually, in the case of the particular derivative, their 
contributions to the vector and tensor potentials Yi and Zij turn out to be zero, 

I^Leibniz^i = , (6.17a) 

^Leibniz^ = , (6.17b) 

while the contribution to the scalar potential T, also in the case of D^'^'^*[F], is found to be 
^ ^ G^m\ , . (\\ , 11G'X"^2 i r,(\ \ , , ^ 

^LeibnizT = ^^>'Ai (^-J + ^ ^12^^ J ^ ^ ^ 

The modification of the acceleration of body 1 which is generated by the latter Leibniz terms 
reads as 



'^Leibniz '^1 
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88 G^niiml 
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12 



(6.19) 
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On the other hand, we computed in Sections |I^ and many distributional terms as- 
sociated with the derivative of the non-linear potentials in the right-hand side of the field 
equations [see ( |3.15| )]. Most of these terms are simply given by the Schwartz distributional 
derivative. The only terms which require the new distributional derivative of |^ come from 
the computation of the "self" parts of the non-compact potentials (see Section |Vp. In this 
case, the modifications of the potentials have been found to be given by ( p.7| ), from which 
we obtain the following modification of the acceleration: 
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Itself a 1 
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(6.20) 
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Adding up ( |6.19| ) and ( |6.20| ) we therefore obtain the total effect of the (particular) distribu- 
tional derivative as 



'^distribution 
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2(ni2t;2)M2 - -vln\2 - ^{ni2V2)v2 



+ 7^^:£F^nl2 . (6.21) 
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Interestingly, this quite simple piece of the acceleration of particle 1 involves the velocity of 
particle 2 alone, and therefore does not stay by itself invariant under a Lorentz transforma- 
tions, or, rather, at this order, a Galilean transformation (indeed, for this to be true the term 
should depend on the relative velocity V12 = vi — V2). Therefore, if we are correct, since we 
are using harmonic coordinates and have employed a Lorentzian regularization, the result 
( |6.21| ) has to combine with other pieces in the acceleration so as to maintain the Lorentz 



invariance of the equations. We have found that this is exactly what happens: the depen- 
dence of ( |6.21| ) over the velocity V2 is mandatory for the Lorentz invariance of the final 3PN 
equations to work. This constitutes, in our opinion, an important check of the relevance of 
the distributional derivative introduced in . It shows also that this derivative is merely a 
tool for preventing a breakdown of the Lorentz invariance when performing integrations by 
parts of complicated divergent integrals [the last term in ( |6.21j ), which is not checked by the 
Lorentz invariance, will turn out to be absorbed into the adjustement of a certain constant, 



see Section VII 



The previous check has been done with the "particular" distributional derivative ( p. 91) , 
and it is interesting to redo the computation in the case of the distributional derivative 
defined by ( p.lO| )- p.ll] ), that we recall is more satisfying than the particular one because 



it obeys the rule of commutation of successive derivatives. (But note in passing that we 
have verified that the particular derivative does not yield any ambiguity at the 3PN order 
which would be due to the non-commutation of derivatives; however, such ambiguities could 
arise at higher orders, in which case the "correct" derivative would be more appropriate.) 
In particular, while the particular derivative is entirely deterministic, the derivative (|2.10|) - 
( p.llp depends on a constant K, and it is important to know the fate of this constant in 
the final equations of motion, and how the test of the Lorentz invariance will manage to be 
satisfied in fine. Like for the particular derivative we find that the incidence of this derivative 
is through two distinct contributions, Leibniz and self. Consider the Leibniz contribution: 
we perform exactly the same computation as before, i.e. based on the formulas in the 
appendix 0, and find that in the case of the new derivative ( p.lO| )-( p.ll| ) the terms in the 
potentials Yi and Zij are no longer zero, but are given by 

^LcibnizK, I = (-^ + G'mlvidij (^^^ +1^2, (6.22a) 
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6 



Leibniz -^ij 



1 

15 



15 



K G-'mid, 



+ 1^2. 



(6.22b) 



Our convention is that the exphcit indication in the left-hand side of the dependence over 
K means that the computation is performed using the "correct" distributional derivative. 
In the case of the modification of the potential T the things are a little more complicated 
because we have to take into account, in addition to a "linear" contribution similar to those 
of ( |ti.22 ), the "non-linear" term which is generated by the modification of the tensor potential 
Zij shown in (|6.22b| ); cf the source term ZydijV in the definition (|3.27a|) of T. We obtain 
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(6.23) 



where u = ^ + ^. Now, using the results of the previous subsection, we see that many 
of these terms are in the form of a gauge transformation corresponding to a gauge vector 
of the type ( |6.11| ). Indeed, we pose 



^.L-4'-->e)'- 
)(■ 



£2 



m 
m J 



With this choice the Leibniz corrections in Yi and Zij become pure gauge, 

'^Leibniz ^ I 
I J 

'^Leibniz -^ij I 



(6.24a) 
(6.24b) 

(6.25a) 
(6.25b) 



while we can re-write ( |6.23| ) in the simplified form 
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(6.26) 



The corresponding modification of the acceleration of particle 1 is found to be 
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(6.27) 
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where the last term represents the gauge term ( |6.15| ) but computed with ( |6.24|) . Therefore, 
modulo a change of gauge, we see that the Leibniz modification of the acceleration brought 
about by the correct derivative has exactly the same form as that, given by ( |6.19| ), due to 
the particular one. However, we must also include the contribution of the self terms. We 
have redone the computation of the self terms as in Section |V| but using the /^-dependent 
derivative and compared the corresponding acceleration with the previous result (|6.2CI|) . We 
get 
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(6.28) 
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Subtracting (|6.19| ) and (|6.27| ) for the Leibniz terms, and adding up the difference of self 
terms given by (|6.28|) , we thereby obtain the difference between the total effects of the two 
distributional derivatives in the acceleration as 



'^distribution'^l I '^distribution'^1 
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(6.29) 
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As we see, there is a dependence on the individual velocity V2 which is left out. Anticipating 
the result that a\, computed with the particular derivative, is invariant under Lorentz trans- 
formations, this means that the ii'-dependent derivative breaks down the Lorentz invariance 
for general values of K (indeed the gauge term cannot modify the behaviour under Lorentz 
transformations). Fortunately, we are now able to fine tune the constant K so that the 
velocity-dependent terms in ( |6.29| ) vanish. Therefore, we obtain a unique value. 



41 

K= , (6.30) 

160 ' ^ ^ 

for which the equations of motion computed with the help of the correct derivative (|2.10|) - 
( p.llp are Lorentz invariant, as they are with the particular derivative. 



Thus, in the case of the correct derivative, the equations of motion are not in general 
Lorentz- invariant, despite the use of the Lorentzian regularization. The likely reason is that 
the distributional derivatives were not defined in a "Lorentzian" way (their distributional 
terms involve the delta-pseudo-function PiSi and not the Lorentzian one Pf Ai). Recall that 
the Lorentzian regularization permitted to add some crucial contributions, proportional to 
m\m2 in the acceleration of particle 1 [see for instance (|5.35|) ], which are mandatory for 



satisfying the Lorentz invariance. In the case of the correct derivative, we find that there is 
still a limited class of terms, proportional to m|, which do not obey the Lorentz invariance, 
unless K is adjusted to the unique value ( |6.30D . Finally we obtain 
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"distribution I — Odistribution^i — 7- ^— E ^12 + Offlx I • (^O.oij 

1^ 10 &rl2 I 

160 160 



We shall see in Section |VII| that the effect of the first term is simply to modify a logarithmic 
constant hi{r'2/ S2) that we shall adjust when we look for a conserved 3PN energy. After 
adjustement of this constant we find that the 3PN equations of motion computed with the 
two derivatives are physically the same since they merely differ by the gauge transformation 
appearing in ( |6.31] ). 



C. Non-distributivity contributions 

The distributive parts of the linear momentum P{ and force Fl densities have been 
written down in ( p.35|) . They were obtained under the uncorrect hypothesis of distributivity, 
that is [-FG]^ = [Fj-j^fG]^, and we must now correct for this. (As explained in Section p| , 
our strategy has been to delineate as much as possible the problems, by concentrating our 
attention first on the computation of the regularized values of the potentials when taken 
individually, second on the corrections due to the non-distributivity, i.e. [FG\ 7^ 
Again, we find that such a subtlety as the non-distributivity makes a difference starting 
precisely at the 3PN order. We get for the required corrections in P{ and Fl: 



+ ^( l2v\[V,V,], - l2v\[V,],\V,],+2v\[V\ 



-?,v\[V],[V\+v\[V]\ 
-S[V^V.],+A[V\m,+A[V]lm, 

+ - \^H[V]\ - lQvi[V^^,], + IQvmAVj]^ , (6.32a) 

+ ^J-vl[VmV],+v![Vd,V],-8[V,Ud,Vj],+8[VAV,l 



+ [V]l[d,V], - [V\[d,V], + 2[V'diV],-2[V],[Vd,Vl 
+ ^e{-r"lVU^^V], + \vt[V^.V], 

+lvl[V]l[d.V], - ^vl[V\[d.Vl+3vl[V'd.Vl-3vl[VUVd,Vl 
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+l[v]l[d.v], - 3[VUV\[d,Vl + -[V\ [d, V] , 
-3[V]l[Vd,V],+2[V\[Vd,V], + 2[Vl[V%Vl - ^[V%Vl 

+8[V],mUVj\, + 8[Vl[d,V,V,], - iQmv,v,i 

+^d.Vl[V,l[Vj], + 4[d,Vl[V,V,], - 8[d.VV,V,l 
-UvlmilVj], + UvfmV.l 
+Avi[V]l[d,Vj], + Avi[V\mi - 8vi[V%V,l 
+8vi[Vl[d,Vl[V,l+8vi[Vd,Vl[V,l - 16vi[Vd,VV,l 
+8vi[Vk],[diWjkl - 8vi[VAW,kl + 8vi[d,Vk],[W,kl - 8vi[d,VkW,kl 
-4:viVl[V]A^^W,kl + 4viv'l[Vd,W,k]^ 
-Avivlld.VllW.kl + 4viv'l[d^VW,k]^ 
+ 16viv'[[d,Vkl[Vjl - lQviv'l[d,VkV,], 

+8[X],[d,V],-8[Xd,V],+8[V],[d,X],-8[Vd,X],] . (6.32b) 



These formulas look complicated but are in fact rather simple to evaluate because they 
require only some lower-order post-Newtonian precision in the potentials, with notably all 
the difficult non-compact potentials needed at the Newtonian order only (hence the interest 
of separating out the problems as we did). Note that it is crucial here to employ the 
Lorentzian regularization [F]-^. The net result of this computation is 
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(6.33a) 
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(6.33b) 



Therefore the supplement of acceleration linked to the non-distributivity is 
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Since only the relative velocity V12 is involved this part of the acceleration is Galilean 
invariant. Furthermore, it can be expressed in a simpler way by introducing an infinitesimal 
gauge transformation of the type (|6.11|) . We pose 
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(6.36) 



Thus, the only Physics brought about by the non-distributivity (i.e. which is not affected by 
a gauge transformation) is constituted by the quartic (G^) term displayed in the right-hand 
side of ( |6.36| ). 



VII. THE 3PN EQUATIONS OF MOTION 



A. Existence of the conserved energy 



At present, the equations of motion are complete. We want now to look for the conserved 
energy associated with these equations at the 3PN order (considering of course only the 
conservative part of the equations, i.e. excluding the radiation reaction acceleration at the 
2.5PN order). We shall see that the existence of an energy is not immediate, but requires 
the adjustment of a certain constant. 

We proved in Section that the equations of motion of body 1 depend on two arbitrary 
constants, which are the constant r'^, tending to zero as we approach the particle 1 (but 
considered here as taking some finite non-zero value), and the constant S2 associated with 
the Hadamard regularization near the other particle 2 [see ( p.3| )]. Similarly, the equations of 
body 2 depend on the constants r'2 and si. All these constants appear inside the logarithms 
entering the equations of motion in harmonic coordinates. Gathering the results for the 
"logarithmic" part of the equations of body 1, we obtain 
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(7.1) 



where the dots indicate the terms which do not contain any logarithms. The terms shown 



in ( |7. Ij ) contain the whole dependence of the acceleration of 1 over r'^ and S2; there are no 
other constants elsewhere. Notice that S2 enters a single quartic-order term proportional to 
G^mimg. Now, most of the terms in ( [7. 1| ) can in fact be gauged away. To see this, we apply 
the formula (|6.15D with the particular choice 
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The corresponding transformation of the acceleration is ( |6.15|) , except that (ei)in and (e2)in 
depend on time through the orbital separation ri2, so in fact this formula should contain 
also some terms proportional to the time-derivatives of (ei)in and (e2)in; but the point for us 
is that these extra terms are free of any logarithms. Therefore, modulo the dots indicating 
the logarithmic-free terms, we can write 



a\ = {S^a\)i^ - -— In (^-^ j + . . . , (7.3) 

where {S^a\)\n denotes the coordinate change of the acceleration. 

The term in ( |7.3| ) which is left out after this coordinate change depends only on the ratio 
between and S2 (similarly, in the equations of motion of body 2, we would find the ratio 
of r[ and Si). This term is of the same type as the one in (|6.31|) giving the difference of 
accelerations, modulo a change of gauge, when different distributional derivatives are used. 
Notice that the constant was originally absent from the equations of motion of 1, but has 
to be introduced in order to "remove" these logarithms by the coordinate transformation. 
Therefore, the only physical freedom remaining in the equations of motion is the yet un- 
specified constant In j . Now we use this freedom to find a conserved energy associated 
with the equations of motion, which means a local-in-time functional E of the trajectories 
and velocities of the two particles, i.e. 



E = i?[yi(t),y2(t);vi(t),V2(t)], (7.4) 
which is constant as a consequence of the 3PN equations of motion, i.e. 



dE ,dE ,dE ,dE .BE —, ^ , ^ 

The accelerations ai and ai are to be replaced by the functionals of the positions and 
velocities given by the 3PN equations of motion. Our special notation for the remainder 
means a radiation-reaction term which is purely of order 2.5PN plus the neglected terms at 
3.5PN; schematically 0(7) = ^Fg + 0(7). See (|7l8|) below for the expression of the term 
-^F^. If an energy exists, the quantity mia\v\ + m2a2f2 must be a total time derivative. In 
practive, we look for a local-in-time functional -D[yi, y2; vi, V2] such that 



mid^v\ + m24^ + = C(7) , (7.6) 

and we obtain the energy as E = |mi + |m2v| + D. Now, the computation with our 3PN 
equations of motion (obtained by means of, say, the particular derivative) shows that the 
quantity D does not exist for any values of the constants In (^2/52) and In (r^/si). However, 
we find that this "nearly" works, because we can determine some ID such that 
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From the computation we obtain D as a well-defined local functional of the positions and ve- 
locities of the particules [containing in particular some logarithms In (ri2/r'^) and In (?"i2/''"2)]- 
The right-hand side of ( [7. 71 ) cannot be written, for generic values of In (^3/52) and In (r'^/si), 
in the form of a total time-derivative. It would be possible, for this to be the case, to adopt 
the simplest choice that both these constants are numerically equal to |||. However, this 
choice does not represent the most general solution for obtaining a total time-derivative. 
Indeed, nothing prevents In {r[/ si) and In (r2/s2) to depend also on the masses mi and m2, 
and therefore such a dependence on the masses should in fact be mandatory (totalitarian 
principle). Since the regularization procedure followed in this paper is more mathematical 
than physical, we can be confident that no Physics will be overlooked only if at each step we 
obtain the most general solution allowed by the process. Unfortunately, the most general 
solution in this case contains an arbitrary parameter. 

The necessary and sufficient condition for the right side of (|7. 7| ) to be a total time- 
derivative is that the factor of {ni2Vi) in (|7.7] ) be invariant by exchanging the particle's 
labels 2 and 1, i.e. 
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Denoting by Am the common value of both sides of 
mi + m2, we obtain the most general solution as 
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81), where A is a constant and m 
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This A is a dimensionless quantity which is the same for the two particles 1 and 2. We now 
prove that A is necessarily a pure nwmencaZ constant, independent of the masses. Notice that 
the A-term in (|7.9a|) will yield a contribution to the acceleration of 1 which is, as concerns the 
dependence over the masses, of the type mim|mA [see (fOp j. If A depends on the masses, 
it must be a symmetric function of mi and m2, and therefore it can be expressed solely in 
terms of the symmetric mass ratio v = "^^^ ■ Suppose that A = Xl-^ AjZ/*, where the Aj's 
are numerical constants, so the A-term in the acceleration of 1 is of the type mim2m ^ AjZ/*. 
First, we see that all the cases i < —I are excluded because the equations of motion would 
not have the correct perturbative limit when z/ — 0; for instance, in the case i = —1, we get 
a term of the type rn?m2 which tends to m2 in this limit, and therefore modifies the geodesic 
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motion of a test particle around a Schwarzschild black hole, which is of course excluded. 
Second, the cases i > I, though they pass the simplest physical requirements, imply that 
the individual particle accelerations are no longer polynomials in the two individual masses 
nil and m2, because of the appearance of inverse powers of the total mass m = mi +7712- For 

2 3 

instance, the case i = 1 leads to a term of the type — But we know that when doing a 
diagrammatic expansion of the N body problem based on the post-Minkowskian expansion 
pOj for the details of the method) that each successive diagram is a polynomial of the 



see 



masses. Therefore, we exclude the possibility that some inverse powers of the total mass 
appear, and find, in conclusion, that A is a pure constant (A = Aq). 

At last, we have succeeded in finding a conserved energy at the 3PN order by specifying 
an unknown logarithmic ratio, but at the price of having introduced an arbitrary purely 
numerical constant A. The constant A will be left undetermined in the present work. So the 
final 3PN equations of motion we obtain in this paper, as well as the final 3PN energy, depend 
on the unknown parameter A. The appearance of A suggests that the present formalism, 
based on a point-mass regularization, is physically incomplete. The resulting ambiguity is 
equivalent to the "static" ambiguity found by Jaranowski and Schafer [Q. It is probably 
linked to the fact that one can write the Einstein field equations into many different forms, 
which are all equivalent in the case of regular sources, but which are in general not equivalent 
in the case of point-particles because the distributional derivative does not obey the Leibniz 
rule. If we had chosen initially a different form of the field equations, the Leibniz terms we 
computed in Section |V|, could have been different. More precisely, only that part of the 
Leibniz terms which is Galilean-invariant and consequently is not required by the Lorentz- 
invariance symmetry could change. But we have seen in ( |6.19| ) and ( |ti.31 ) that the Galilean- 
invariant part of the Leibniz terms is precisely made of only one term, which is of the same 
type (proportional to G^mim2) as the term containing the constant In(r2/s2) [see ( [7.3|) ] that 
we have adjusted in the equations (|7.9a|) , resulting in the appearance of the constant A. Thus, 
in agreement with Jaranowski and Schafer |Q, we might say that A encodes an ambiguity 
associated with the violation of the Leibniz rule by the distributional derivative. At a deeper 
level, this would mean that the ambiguity is a consequence of a theorem of Schwartz ||6^ 



according to which it is impossible to define at once a multiplication of distributions which 
agrees with the ordinary product for continuous functions, and a derivation of distributions 
which satisfies the Leibniz rule and reduces to the ordinary derivative in the case of 
functions. If this explanation is correct, it is unlikely that the constant A could be determined 
within the present formalism. 

We find by combining (|7.9|) and (|7.7|) that the dependence of the 3PN energy £" on A is 



E = E-'-lx?^4^, (7,10) 

where E does not depend on A, while we obtain, using ( [7.3|) , that the acceleration writes 
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44 , G^mim'^m 

a A ^ r n 

^ 3 



12 ' 



(7.11) 
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where similarly d\ is independent of A. On the other hand, the acceleration and energy 
depend also on the two constants r[ and r'^, but from the previous discussion this is not a 
problem because r[ and are associated with an arbitrariness in the choice of coordinates: 
the 3PN equations of motion contain the logarithms ln{ri2/r[) and In(ri2/r2) which have 
been shown in (|7.1|) to be in the form of the gauge transformation associated with ( [7.21) . 
In particular, the "constants" Inr'^ and Inrg, which might be said to be formally infinite 
because r[ and r'2 were tending to zero [recall the discussion after ( p.l^ )], will never appear 
in any physical result. Similarly, the dependence of the energy on the logarithms In {ri2/r[) 
and In (ri2/r2) is pure gauge. From (|6.16|) we get 
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+ 3(ni2fl)(ni2^^12) - (^^1^^12) 



ln(^ 



+ 1 



(7.12) 



where the dots denote the terms independent of logarithms [this result can also be checked 
directly using ([7.11)]. 

Finally, to be more specific about the influence of the distributional derivative, notice 
that the solution we have obtained in ( [7.9[) corresponds to the "particular" distributional 
derivative defined by ([2.9| ). If one uses the "correct" derivative ( p.lO|) -( prTT| ) instead, with 
the value K = we have obtained in (|6.3CI|) from the Lorentz invariance, we obtain the 
same equation to be solved as ( [7.7[) but with the rational fraction 
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instead of 



This is easily seen thanks to (|6.31[) . So, the solution becomes in this case 
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(7.13a) 
(7.13b) 



Replacing this into the equations of motion (and associated energy), it is then clear that they 
are physically the same as those computed with the other derivative, because they differ by 
the mere change of gauge. 



41 
160 



41 
160 



that we obtained in ( [6.31[ ). We give it here thoroughtly for completeness: 
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(7.14) 



(7.15) 
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B. End results 



We present the 3PN equations of motion of the particle 1 in harmonic coordinates, which 
are obtained by summing up all the contributions of the potentials computed in Sections ^ 
and [V|, as well as the pieces due to the non-distributivity and the Leibniz terms (see Section 
|VID . The equations depend on two gauge constants r'^ and r'2 through some logarithms, and 
on one unknown purely numerical coefficient A. The equations of the particle 2 are obtained 
by exchanging all the labels 1^2. 
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These equations are in full agreement with the known results valid up to the 2.5PN or- 
der ||2^- p6i^ . They have the correct perturbative limit given by the geodesies of the 
Schwarzschild metric at the 3PN order. Most importantly, the equations are invariant un- 
der Lorentz transformations (developed to 3PN order); this can be checked using for instance 
the formulas developed in |^8|. Finally, as we have seen previously, the equations of motion 
admit a conserved energy at the 3PN order. The study of the Lagrangian (and Hamiltonian) 
formulation of these equations is reported in a separate work |37|. The energy is given by 
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This energy is conserved in the sense that its time- derivative computed with the 3PN equa- 
tions of motion equals the radiation reaction effect at the 2.5PN order, namely 
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The rather complicated expressions (|7.16|) - (|7.17|) simplify drastically in the case where the 
orbit is circular [apart from the gradual inspiral associated with the balance equation ( [7.18|) ] 
and where we place ourselves in the center-of-mass frame. The circular orbit corresponds 
to the physical situation of the inspiralling compact binaries which motivate our work pl- 



Here, we give the result concerning circular orbits without proof (see also |0). The relative 
acceleration reads 



dwi2 , .2 
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wVl2 + Freac + 0(7) , (7.19) 
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where Fj-eac is the standard radiation-reaction force in harmonic coordinates, 
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being the symmetric mass ratio), and where the orbital frequency uj of the relative 



circular motion is given to the 3PN order by 
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7, and we recall that ri2 = |yi — y2| is 



The post-Newtonian parameter is defined by 7 
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the orbital separation in harmonic coordinates. The constant Tq appearing in the logarithm 
is related to the two constants r[ and by 
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The 3PN energy E in the center of mass of the particles, which is such that ^ = as a 
consequence of the conservative equations of motion, is obtained as 
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The invariant 3PN energy follows from the replacement of the post-Newtonian parameter 
7 by its expression in terms of the frequency uj as deduced from computing the inverse of 
{ \ni\ ). We find 



E = ——fic^x-l^ 1 



3 1 \ / 27 19 1 
U \ X + 1 u 

4 12 / V 8 8 24 



2 1 2 

U X 



+ 



675 



64 



+ 



209323 205 o HO, 



4032 



96 



155 
"96"' 



35 
5184" 



iy']x-' + 0{8) 



(7.24) 



where the parameter x is defined by 



(7.25) 
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Note that the logarithm disappeared from the invariant expression of the energy ( [7.24|) , in 
agreement with the fact that it is pure gauge. However, the constant A stays in the final 
formula; the static ambiguity constant uOg of Jaranowski and Schafer is related to it by 



11 \ _ 1987 
3 840 



see 
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APPENDIX A: SUM OF LEIBNIZ TERMS 

In this appendix we give the sum of all the terms of the type ^Leibniz^ introduced by 
( p. 25 ) that we have encountered during the process of simplification of the 3PN potentials. 



The reduction of these terms using the distributional derivative is done in Section W\ 
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